Control Moment Gyroscope Skew Angle Variation and Singularity Penetration

Abstract

This manuscript contains a brief introduction of Control Moment Gyroscopes followed by the mathematical basis for potential singularities and an analysis of how skew angle variability may impact their occurrence. MATLAB is utilized as the primary modeling tool along with WolframAlpha for mathematical derivations of matrix determinants. The results of the modeling efforts show that a uniform skew angle of 90° allows a high maximum angular momentum. Additionally, we attempt to show that having two CMGs at a skew angle of zero could result in similar gains as a uniform 90° configuration and briefly introduce singularity penetration.

1. Introduction

The use of non-tumbling spacecraft, of either the spin stabilized or 3-axis stabilized variety, imposes a requirement for proper attitude control. Control Moment Gyroscopes (CMGs) are pervasive in fulfillment of this requirement not only because of the higher torque profile afforded over a reaction wheels but also due to greater power efficiency gained through their use [1]. These moment exchange devices utilize the principle of conservation of momentum to bypass the fuel and weight limitations of thrusters while still being able to provide a fairly significant slew rate for spacecraft attitude control. This level of performance does not come without its own share of drawbacks: momentum saturation, bearing decay, and most notably singularities all impose operational constraints [2].

CMG steering was implemented into a Simulink model as shown in Figure 1 in order to analyze these singularities. This model also contained sinusoidal trajectory generation, system dynamics with feedforward control, and a combination of quaternions and direction cosine matrices for the rotational kinematics. The simulations allowed us to observe that varying CMG geometric configuration could have a significant impact on singularities occurrences. With this knowledge in hand one can reason that a spacecraft designer can and should take into account CMG skew variations in order to optimize maneuvering for particular applications. For instance, a communication satellite in geosynchronous orbit and a low earth orbit imaging satellite could end up with very different attitude, pointing, and slew rate problems to solve.

Figure 1.

Simulink model topology.

A visual representation of CMG placement can be seen in the two depictions of Figure 2, which were adapted from drawings in Ref. [3].

Figure 2.

(A) and (B) CMG arrangement.

Figure 2(A) shows a typical non-redundant skewed pyramid design with a uniform skew angle ( $\beta$ ). While typical space applications maintain constant skew angles amongst all the CMGs, it is possible to implement mixed a slew angle design as demonstrated in Figure 2(B) where CMG #3 has a different skew angle (a roughly 90° skew angle is shown). These figures also depict another critical variable, gimbal angle. Gimbal angle, denoted by $\delta$ in Figure 2 and $\theta$ in the equations of section II, is used to show the rotation that the CMG undertakes in order to impart a torque on the space vehicle. When CMGs are gimbaled, the change in their position alters momentum vector direction (h₁, h₂, and h₃). This in turn changes the torque vector which is perpendicular to both the gimbal axis and the momentum vector. For instance, CMG #3 in depiction (B) will have torque in the shaded plane and it is when a torque is required outside of that plane that a singularity may occur.

These fundamental ideas form the basis for CMG operation and the results shown later in Section III through Section V. We will attempt to show that mixed skew angle CMG configurations may provide some flexibility when dealing with singularities and briefly illustrate some applications and limitations for singularity penetration.

2. Mathematical background

Before going into the results of skew angle variation, let us first set the stage for our future analysis with mathematical derivations used to form the basis for CMG steering logic. For more detailed background information, consult Refs [4,5] for an in depth treatment of kinematics and trajectory generation and Ref. [6] for a complete overview of material found in this section. In order to translate a commanded spacecraft movement into desired torque and CMG gimbal angles rates, we derived the mathematical basis for momentum vectors (h_x, h_y, and h_z) from Figure 2(A) and (B) as shown in Eq. (1) through Eq. (3).

$$h_x=-cos{\theta }_1+\mathit{\text{cosβsin}}{\theta }_2+\mathit{cos}{\theta }_3$$ E1

$$h_y=-\mathit{\text{cosβsin}}{\theta }_1-cos{\theta }_2+\mathit{\text{cosβcos}}{\theta }_3$$ E2

$$h_z=\mathit{\text{cosβsin}}{\theta }_1+\mathit{\text{cosβsin}}{\theta }_2+\mathit{\text{cosβsin}}{\theta }_3$$ E3

These momentum vectors were then used to obtain a matrix representing the spatial gradient of torque ( $\frac{\partial H}{\partial \theta }$ ) by taking the partial derivative of each component with respect to gimbal angle. This derived matrix is referred to as the [A] matrix and is shown below in Eq. (4). The rows correspond to h_x, h_y, and h_z components and the columns are the partial derivative of the momentum components with respect to each gimbal angle.

$$H=\frac{\partial H}{\partial \theta }=\left[\begin{array}{ccc}\mathit{sin}{\theta }_1& \mathit{\text{cosβcos}}{\theta }_2& -\mathit{sin}{\theta }_3\\ -\mathit{\text{cosβcos}}{\theta }_1& \mathit{sin}{\theta }_2& \mathit{\text{cosβcos}}{\theta }_3\\ \mathit{\text{sinβcos}}{\theta }_1& \mathit{\text{sinβcos}}{\theta }_2& \mathit{\text{sinβcos}}{\theta }_3\end{array}\right]$$ E4

Torque was found via a change in angular momentum as shown in Eq. (5) and finally a desired gimbal rate to produce this torque was found via Eq. (6).

$$\dot{H}=\frac{\partial H}{\partial \theta }\frac{\partial \theta }{\partial t}=\left[A\right]\dot{\theta }$$ E5

$${\left[A\right]}^{-1}\dot{H}={\left[A\right]}^{-1}\left[A\right]\dot{\theta }=\dot{\mathrm{\theta }}$$ E6

Upon closer analysis we can see that there is potential for the ${\left[A\right]}^{-1}$ matrix in Eq. (6) to become singular when applying the matrix inverse per Eq. (7) [7].

$${\left[A\right]}^{-1}=\frac{1}{\mathit{det}\left[A\right]}\left[\mathit{\text{Cofactors of}}\left[A\right]\mathit{\text{matrix}}\right]$$ E7

If the determinant of the [A] matrix becomes zero, then we have a situation where $\frac{1}{\mathit{det}\left[A\right]}$ is equal to $\frac{1}{0}$ leading to failures in the CMG steering logic. Section III will dig deeper into this singularity analysis.

3. Singularity analysis

This section reemphasizes and builds on work presented in Refs [8, 9]. To begin to analyze when CMG Singularities occur, we must first look into the factors which drive their occurrence. Specifically, we know that when the determinant of the [A] matrix is equal to zero it means that at least one or a combination of multiple CMGs are incapable of producing the required torque. We turned to the WolframAlpha website to solve for the determinant of the [A] matrix from Eq. (4) and produce Eq. (8) [10].

$$\begin{array}{l}det\left[\mathrm{A}\right]=sin\left(\mathrm{\beta }\right)\left\{2{\mathbf{cos}}^{\mathbf{2}}\left(\mathbf{\beta }\right)\mathbf{cos}{\mathbf{\theta }}_{\mathbf{2}}cos{\mathrm{\theta }}_{1}cos{\mathrm{\theta }}_3-\mathbf{cos}\left(\mathbf{\beta }\right)\mathbf{cos}{\mathbf{\theta }}_{\mathbf{2}}sin{\mathrm{\theta }}_{1}cos{\mathrm{\theta }}_3\right.\\ \left.+\mathbf{cos}\left(\mathbf{\beta }\right)\mathbf{cos}{\mathbf{\theta }}_{\mathbf{2}}cos{\mathrm{\theta }}_{1}sin{\mathrm{\theta }}_3+sin{\mathrm{\theta }}_1\mathbf{sin}{\mathbf{\theta }}_{\mathbf{2}}cos{\mathrm{\theta }}_3+cos{\mathrm{\theta }}_1\mathbf{sin}{\mathbf{\theta }}_{\mathbf{2}}sin{\mathrm{\theta }}_3\right\}\end{array}$$ E8

Combining like terms allows us to factor out a $cos\left(\mathrm{\beta }\right)cos{\theta }_2$ term from the first three lines and apply Eq. (9) for terms with ${\theta }_1$ and ${\theta }_3$ . Likewise we could factor a $sin{\theta }_2$ term out of the last two lines and apply Eq. (10) to ${\theta }_1$ and ${\theta }_3$ .

$$sin\left(\mathrm{A}-\mathrm{B}\right)=\text{cosBsinA}-\text{sinAcosB}$$ E9

$$sin\left(\mathrm{A}+\mathrm{B}\right)=\text{cosBsinA}+\text{sinAcosB}$$ E10

The final result is Eq. (11) shown below where we have abbreviated sine and cosine functions as “s” and “c” respectively.

$$det\left[\mathrm{A}\right]=sin\left(\mathrm{\beta }\right)\left\{\mathrm{s}{\mathrm{\theta }}_2\left[\mathrm{s}\left({\mathrm{\theta }}_1+{\mathrm{\theta }}_3\right)\right]+\mathrm{c}\left(\mathrm{\beta }\right)\mathrm{c}{\mathrm{\theta }}_2\left[\mathrm{s}\left({\mathrm{\theta }}_3-{\mathrm{\theta }}_1\right)+2\mathrm{c}{\mathrm{\theta }}_1\mathrm{s}{\mathrm{\theta }}_3\mathrm{c\beta }\right]\right\}$$ E11

Analyzing Eq. (11) leads us to several cases where singularities can occur as shown in Eqs. (12)–(17). A bolded strike-through represent a condition where the altered terms are equal to zero in order to make the right side of Eq. (11) zero (i.e., singularity conditions).

$$\mathit{s}\left(\mathbf{\beta }\right)\left\{\mathrm{s}{\mathrm{\theta }}_2\right[\mathrm{s}\left({\mathrm{\theta }}_1+{\mathrm{\theta }}_3\right)+\mathrm{c}\left(\mathrm{\beta }\right)\mathrm{c}{\mathrm{\theta }}_2\left[\mathrm{s}\left({\mathrm{\theta }}_3-{\mathrm{\theta }}_1\right)+2\mathrm{c}{\mathrm{\theta }}_1\mathrm{s}{\mathrm{\theta }}_3\mathrm{c\beta }\right]\}$$ E12

$$s\left(\mathrm{\beta }\right)\left\{\mathrm{s}{\mathrm{\theta }}_2\right[\mathbf{s}\left({\mathbf{\theta }}_{\mathbf{1}}+{\mathbf{\theta }}_{\mathbf{3}}\right)+\mathrm{c}\left(\mathrm{\beta }\right)\mathrm{c}{\mathrm{\theta }}_2\left[\mathbf{s}\left({\mathbf{\theta }}_{\mathbf{3}}-{\mathbf{\theta }}_{\mathbf{1}}\right)+\mathbf{2}\mathbf{c}{\mathbf{\theta }}_{\mathbf{1}}\mathbf{s}{\mathbf{\theta }}_{\mathbf{3}}\mathbf{c\beta }\right]\}$$ E13

$$s\left(\mathrm{\beta }\right)\left\{\mathrm{s}{\mathrm{\theta }}_2\right[\mathrm{s}\left({\mathrm{\theta }}_1+{\mathrm{\theta }}_3\right)+\mathrm{c}\left(\mathrm{\beta }\right)\mathrm{c}{\mathrm{\theta }}_2\left[\mathrm{s}\left({\mathrm{\theta }}_3-{\mathrm{\theta }}_1\right)+2\mathrm{c}{\mathrm{\theta }}_1\mathrm{s}{\mathrm{\theta }}_3\mathrm{c\beta }\right]\}$$ E14

$$s\left(\mathrm{\beta }\right)\left\{\mathbf{s}{\mathbf{\theta }}_{\mathbf{2}}\right[\mathrm{s}\left({\mathrm{\theta }}_1+{\mathrm{\theta }}_3\right)+\mathbf{c}\left(\mathbf{\beta }\right)\mathbf{c}{\mathbf{\theta }}_{\mathbf{2}}\left[\mathrm{s}\left({\mathrm{\theta }}_3-{\mathrm{\theta }}_1\right)+2\mathrm{c}{\mathrm{\theta }}_1\mathrm{s}{\mathrm{\theta }}_3\mathrm{c\beta }\right]\}$$ E15

$$s\left(\mathrm{\beta }\right)\left\{\mathrm{s}{\mathrm{\theta }}_2\right[\mathbf{s}\left({\mathbf{\theta }}_{\mathbf{1}}\mathbf{+}{\mathbf{\theta }}_{\mathbf{3}}\right)+\mathbf{c}\left(\mathbf{\beta }\right)\mathbf{c}{\mathbf{\theta }}_{\mathbf{2}}\left[\mathrm{s}\left({\mathrm{\theta }}_3-{\mathrm{\theta }}_1\right)+2\mathrm{c}{\mathrm{\theta }}_1\mathrm{s}{\mathrm{\theta }}_3\mathrm{c\beta }\right]\}$$ E16

$$s\left(\mathrm{\beta }\right)\left\{\mathbf{s}{\mathbf{\theta }}_{\mathbf{2}}\right[\mathrm{s}\left({\mathrm{\theta }}_1+{\mathrm{\theta }}_3\right)+\mathrm{c}\left(\mathrm{\beta }\right)\mathrm{c}{\mathrm{\theta }}_2\left[\mathbf{s}\left({\mathbf{\theta }}_{\mathbf{3}}-{\mathbf{\theta }}_{\mathbf{1}}\right)+\mathbf{2}\mathbf{c}{\mathbf{\theta }}_{\mathbf{1}}\mathbf{s}{\mathbf{\theta }}_{\mathbf{3}}\mathbf{c\beta }\right]\}$$ E17

Note that some of these cases only drive singularities under limited ranges. For instance, Eq. (12) only drives the determinant to zero when sine $\mathrm{\beta }$ is equal to zero; i.e., when $\mathrm{\beta }$ = 0° (i.e., at a single point skew variation). Additionally, Eq. (16) is only zero when either cosine $\mathrm{\beta }$ or cosine ${\mathrm{\theta }}_2$ is equal to zero; i.e., when either $\mathrm{\beta }$ or ${\mathrm{\theta }}_2$ is equal to 90°.

When the collective effects of Eqs. (12)–(17) are taken into account we obtain Figure 3 which shows us the maximum magnitude of angular momentum that we can achieve for a given skew angle without running into a singularity. This plot was formed from a discretization of roughly 1 million points for curve smoothing purposes.

Figure 3.

Maximum angular momentum achieved before encountering singularity (uniform skew angle).

Here we found a maximum allowable Momentum which corresponds to β = 90°. At this point you can achieve one CMG worth of angular momentum without a possibility of singularity occurrence.

Note that Figure 3 corresponds to Figure 2(A) where all three CMGs have the same skew angle. This graph indicates that when a uniform skew angle of 90° is used, called a box or roof configuration, the spacecraft can exert the highest torque values without experiencing a singularity. Similar results were found in [8] where the researchers indicated that this uniform skew angle of 90° could provide “far superior performance than the commonplace 4 CMG pyramid skewed at 54.73”.

Note that it is also possible to operate in a condition similar to Figure 2(B) where the CMGs can have non-uniform beta angles. The results of mixed skew angle testing are shown in Figure 4 with its conditions tabulated in Table 1.

Figure 4.

Singularity surfaces by skew angle.

	CMG#1 skew	CMG#2 skew	CMG#3 skew	Preferred direction
Case 1	0°	90°	90°	N/A
Case 2	90°	0°	90°	Z
Case 3	90°	90°	0°	N/A
Case 4	0°	0°	90°	X
Case 5	0°	90°	0°	Y
Case 6	90°	0°	0°	X

Table 1.

Mixed skew angle test cases and preferred direction.

*Skew angle β is shown in degrees.

From Figure 4, you can see that some mixed skew conditions lend themselves to certain applications. For instance, Case 2 allows a similar amount of maximum angular momentum in the X and Y directions but a much greater amount in the Z direction. Both Case 4 and 6 allow more movement in the X direction while Case 5 allows more freedom in the Y direction. Interestingly, note that Case 1 and 3 appear to have no clear advantage like the other four cases. In fact, the only case with two 90° skew angles that has a clear preferred direction is Case 2 while the cases with two 0° skew angles always seem to have a preferred direction. This conclusion is backed up by further analysis conducted in Section IV.

4. Further mixed skew insight

In Section III, we saw that some cases of mixed skew angle CMG configurations could have drastic benefits when moving in a particular direction. This begs the question, what is the maximum angular momentum value that can be achieved in a mixed skew configuration and is there a specific configuration that solves all singularity problems. In order to investigate this we attempted to recreate Figure 3 but instead of altering a single, uniform skew angle like Figure 2(A) we iterated through possible skew angle combinations by increments of five. This involved running a set of about 7000 combinations of skew angle through our MATLAB model. The model analyzed the impact of skew angles on the maximum angular momentum that could be achieved without encountering a singularity. The results of this are plotted in Figure 5.

Figure 5.

Maximum angular momentum achieved before encountering singularity (mixed skew angle).

Note that once again the maximum amplitude is one CMG worth of angular momentum but also that there is an almost noticeable pattern to the results. In the results of Figure 5, we can see that iterations 201, 301, 401, 501, 601, 701, and 801 all cause spikes in the maximum allowable torque. The skew angle combinations responsible for these spikes can be seen in Table 2.

Iteration	CMG#1 skew	CMG#2 skew	CMG#3 skew
201	20	0	0
301	30	0	0
401	40	0	0
501	50	0	0
601	60	0	0
701	70	0	0
801	80	0	0

Table 2.

Mixed skew angle test cases and preferred direction.

These results indicate that having two CMGs with a skew angle of zero results in high allowable angular momentum values. It should also be noted that in Figure 5 there are smaller peaks that almost reach an angular momentum value of 1. For example, Figure 5 iteration 704 corresponds to a CMG configuration where CMG#1 is 5°, CMG#2 is 90°, CMG#3 is 0°, allowing a maximum angular momentum of about 0.9. Space is at a premium on a spacecraft and this information could prove useful if certain CMG configuration can be made to fit onboard.

5. Singularity penetration

While the singularities themselves do prove to be a limiting factor it must be noted that there are various techniques available to mitigate their effects. For instance, Refs [6, 11] detail methods of singularity avoidance. A newer technique, illustrated in Ref. [12] espouses a system of singularity penetration called Singularity Penetration with Unit Delay (SPUD). The topology of a circuit for this implementation is depicted in Figure 6.

Figure 6.

Singularity penetration with unit delay circuit.

This circuit takes in the $\dot{\mathrm{\theta }}$ term as seen in Eqs. (5) and (6). These were implemented within the CMG Steering block of Figure 1’s topology. In the topology of Figure 6, a singular condition is monitored and the input from the Inverse Condition Number block is sent in as a threshold value of the switching circuit. Once a specified threshold is reached the circuit will hold onto the $\dot{\mathrm{\theta }}$ value until the singular condition passes. A truly singular condition (e.g., dividing by zero) usually only occurs during one time step, but that one instant is enough to induce extreme errors into a system. We conducted testing under a mildly singular condition where the inverse condition number was only on the order of 1e-6 with the SPUD circuit engaging when it detected the inverse condition number was below 1e-4. The results were not compelling enough to illustrate here due to not being sufficiently singular. However, when a highly singular condition is encountered, SPUD can demonstrate drastic improvements. The results of Ref. [12], tested under a highly singular condition, show a decrease in average Euler Angle tracking error as shown in Table 3.

	Mean error magnitude
System state	X-axis	Y-axis	Z-axis
SPUD inactive	4.4321	5.1939	4.3712
SPUD active	0.0674	0.1063	0.4682

Table 3.

Euler tracking error reduction.

While SPUD provides sizable gains it should be noted that there are factors that must be taken into account when it is utilized. The singular conditions shown in Eq. (12), are triggered when $\mathrm{\beta }$ is equal to 90°. The circuit of Figure 1 was implemented with a step size of 0.0001 and a Runge-Kuta Solver with a variable $\mathrm{\beta }$ input and a yaw maneuver of around 60°. The $\mathrm{\beta }$ was implemented with a slope circuit that ran from −90° to +90°. When the slope circuit crossed the 0° mark a minimum inverse condition number on the order of 1e-20 was reached; well below machine epsilon. The results of this simulation are shown in Figure 7. SPUD was activated at a threshold of 5e-6 in both cases but the maneuver time was implemented as 10 seconds for (A) and 20 seconds for (B).

Figure 7.

(A) and (B) SPUD results with variable $\beta .$

In depiction (A) SPUD reduced average error by about 1e-3 but in depiction (B) SPUD actually induced more error in the yaw channel (i.e., the maneuvering axis) by over five times. Results where SPUD reduced error in the auxiliary, non-maneuvering roll and pitch axes while increasing error in the maneuver axis of yaw were common when SPUD was utilized with a sliding $\mathrm{\beta }$ angle. The sliding $\mathrm{\beta }$ input forces the SPUD circuit into a prolonged singular condition under which it is less effective. The reason for this is that SPUD is designed to penetrate a singular point encountered in a guidance system instead of operate around a singular point itself. Tools are effective when utilized for their specific purpose and usually less effective when implemented improperly; SPUD is no exception.

6. Future efforts

While it was found that having two CMGs with the same skew angle could help achieve a higher allowable torque, there is still much work to be done. Future efforts will look at a deeper resolution (skew angle increments less than 5 degrees) and try to make mathematical sense of why the results are as seen in Figure 5. Additionally, while singularity penetration has proven to be a highly useful alternative to singularity avoidance, further testing should be conducted to analyze when it should be activated. The ideal value of inverse condition number under which SPUD should be activated or ideal rate of change of the condition number are primed for investigation. The key is finding the correct way to identify precursor conditions for an acute singularity occurrence.