Modeling and Analysis of Micro-Spacecraft Attitude Sensing with Gyrowheel

Abstract

This paper proposes two kinds of approaches of angular rate sensing for micro-spacecraft with a gyrowheel (GW), which can combine attitude sensing with attitude control into one single device to achieve a compact micro-spacecraft design. In this implementation, during the three-dimensional attitude control torques being produced, two-dimensional spacecraft angular rates can be sensed from the signals of the GW sensors, such as the currents of the torque coils, the tilt angles of the rotor, the motor rotation, etc. This paper focuses on the problems of the angular rate sensing with the GW at large tilt angles of the rotor. For this purpose, a novel real-time linearization approach based on Lyapunov’s linearization theory is proposed, and a GW linearized measurement model at arbitrary tilt angles of the rotor is derived. Furthermore, by representing the two-dimensional rotor tilt angles and tilt control torques as complex quantities and separating the twice periodic terms about the motor spin speed, the linearized measurement model at smaller tilt angles of the rotor is given and simplified. According to the respective characteristics, the application schemes of the two measurement models are analyzed from the engineering perspective. Finally, the simulation results are presented to demonstrate the effectiveness of the proposed strategy.

1. Introduction

Presently, more and more researchers are focusing on the related technologies of micro-spacecraft [1,2,3]. For micro-spacecraft, the attitude control system (ACS) is one of the major contributors to the mass, volume, power and cost [4]. It will be of great significance for the development of micro-spacecraft if the realization of spacecraft angular rate sensing and the output of three-dimensional control torques can be simultaneously conducted in one instrument [5,6]. For this purpose, some innovative actuators like a variable speed control moment gyroscope (VSCMG) [7], a tilted wheel [8,9], and so on [10], which also have potential to realize the function of the sensors in principle, have been developed. Due to the existence of nonlinear friction from the support structure of some actuators, it is complicated to realize the angular rate sensing for them. However, besides that, according to the difference of the gyro rotor support and combining the advantages of the actuator-CMG [11] with the ones of the sensor dynamically-tuned gyroscope (DTG) [12], the integrated devices can be roughly divided into two categories: the magnetically-suspended double-gimbal control moment gyroscope (MSDGCMG) [13] and the integrated device-based flexible gimbal support structure represented by the gyrowheel (GW) [14]. The former one supports the rotor by active magnetic bearings (AMBs). The latter, GW, is developed based on the principle of DTG by Bristol Aerospace Company for the Canadian Space Agency’s SCISAT-1 Scientific in 2003 [15], and the GW rotor is supported by crossed torsion springs and a gimbal. The MSDGCMG implements the functions of two degrees of freedom (DOF) torque output and two-axis angular rate sensing through working in two different operation modes [16]. However, due to the complicated structure, high mass and large volume, the application of MSDGCMG in micro-spacecraft is seriously restricted. Unlike MSDGCMG, the GW not only can achieve three-dimensional torque output and two-axis angular rate sensing at the same time, but also has the advantage of being light weight and low power, which brings the hope that the GW will have promising application prospects. Therefore, it makes sense to develop the method of the angular rate sensing with GW as an actuator and a redundant measurement device simultaneously.

GW and DTG have similar structures. However, DTG always operates in the fixed tuned speed, and the spin axis of the rotor cannot tilt in radial directions theoretically. However, GW is almost one thousand times greater than DTG in mass, moment of inertia and angular momentum [17]. To realize torque outputs along the radial directions of the rotor, the tilt range of the spin axis of the rotor in GW is significantly increased up to 7° [18]. To realize torque output along the spin direction, the operating speed of the rotor is always kept as time varying. Obviously, GW has more complex dynamical characteristics than DTG. Thus, it is more complicated for GW to realize the spacecraft angular rate sensing while three-dimensional control torques are outputted.

In order to realize the angular rate sensing of the carrier by GW, Dr. Own at Carleton University in Canada realized his work by linearizing the equations of motion of GW [19] at zero tilt angles of the rotor. However, a higher measurement accuracy can be obtained only if the tilt angles of the rotor spin axis are limited to a very small range. This disadvantage imposes strict restrictions on the capacity of the GW output torque, while a higher measurement accuracy is needed. Although Jeffrey M. Hall at Carleton University improved the measurement accuracy of the two-dimensional rate sensing through the ground calibration for the GW [20], the accuracy loss caused by linearization at zero tilt angles cannot be easily compensated at larger tilt angles. Moreover, before the linearization of the equations of motion of GW, an assumption that motor spin speed was always kept constant was made by Dr. Own, which means that the spacecraft angular rates cannot be accurately measured while the control torques along the motor spin axis are outputted by adjusting the motor spin speed. Liu proposed the nonlinear algebraic measurement method of the spacecraft angular rates with the GW supported by the torsion springs through analyzing the full GW dynamics equations [21], but it is hard to calibrate the compensable gyroscopic drift for this nonlinear algebraic measurement method [22], which directly leads to the inaccuracy of the measurement results. Liu also proposed a dynamic measurement method by establishing GW nonlinear state equations [23]. However, it is also difficult to compensate the modeling errors due to the nonlinear form of the state equations.

To overcome the aforementioned drawbacks of the spacecraft angular rate sensing with GW at present, the reminder of this paper is divided into five sections: In Section 2, the GW is succinctly described, and then, its dynamics equations expressed by case coordinates are developed by the Lagrange equations of the second kind. In Section 3, the GW real-time Lyapunov linearization measurement model at arbitrary operating points is established based on Lyapunov’s linearization theory. Especially, by the complex quantity method, the GW complex differential equations within small tilt angles of the nominal position are derived. The small tilt measurement model is established and simplified by ignoring the twice periodic components about the motor spin rate. In Section 4, according to the characteristics of these two proposed measurement models, the application schemes are analyzed from the perspectives of the static measurement and dynamic measurement. In Section 5, simulations are performed to illustrate the validity of the proposed linearization measurement models and their application schemes in this paper. Finally, we draw the conclusions in Section 6.

2. Descriptions of Gyrowheel

2.1. Gyrowheel Physical Structure

The computer-aided design diagram and simplified structure of the GW are shown in Figure 1a,b, respectively. The GW system mainly consists of the following subassemblies: case, motor, flexible suspension structure, flywheel rotor, torquer consisting of current coil and permanent magnet and tilt sensor. Among them, the case is fixed on the carrier, such as spacecraft. The flexible suspension structure is made up of gimbal, inner and outer torsion springs, as shown in Figure 1b; the gimbal is connected to the motor shaft by a pair of inner torsion springs, and the rotor is connected with the gimbal by a pair of outer torsion springs. The rotor driven by the brushless DC motor rotates in a high time-varying speed. Thus, the torque along the spin direction of the rotor can be generated by adjusting the motor speed. Two pairs of torquers perpendicular to each other can provide two-dimensional tilt control torques to make the spin axis of the rotor tilt along the radial directions. Due to the existence of the angular momentum, the radial control torques can be outputted based on the CMG principle. The tilt sensors are designed to measure the tilt angles of the rotor relative to the case. The special physical structure of GW determines that the device can realize the functions of the sensor, like DTG, and the actuator, like VSCMG, at the same time.

2.2. Equations of Motion for an Idealized Gyrowheel

Four body frames and three generalized coordinates are given in Figure 1b for deriving GW dynamics equations using Lagrange’s method. The four body frames are the case frame (${\mathcal{F}}_0$:O-$x_c$$y_c$$z_c$), the motor body frame (${\mathcal{F}}_1$:O-$x_m$$y_m$$z_m$), the gimbal body frame (${\mathcal{F}}_2$:O-$x_g$$y_g$$z_g$) and the rotor body frame (${\mathcal{F}}_3$:O-$x_r$$y_r$$z_r$), respectively. The coordinates (${\theta }_x,{\theta }_y,{\theta }_z$) are the defined generalized coordinates, where ${\theta }_x$ and ${\theta }_y$ represent the rotation angles of the inner and outer torsion springs, respectively, and ${\theta }_z$ represents the spin angle of the motor shaft. The relationship between these four body frames and the generalized coordinates can be further expressed by Figure 2.

Assume that the angular velocity of the GW case fixed on the spacecraft, with respect to inertial space in the case frame ${\mathcal{F}}_0$, is ${\mathit{\omega }}_{\mathit{b}}$ = ${\left[{\omega }_{bx}{\omega }_{by}{\omega }_{bz}\right]}^T$.

Based on the relationship in Figure 2, the angular velocity of the motor shaft ${\mathit{\omega }}_{\mathit{m}}$ in the motor body frame ${\mathcal{F}}_1$ is described as the following equation:

$${\mathit{\omega }}_{\mathit{m}}=\left[\begin{array}{c}{\omega }_{mx}\\ {\omega }_{my}\\ {\omega }_{mz}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ {\dot{\theta }}_z\end{array}\right]+{({\theta }_z)}_z·{\mathit{\omega }}_{\mathit{b}}=\left[\begin{array}{c}{\omega }_{bx}{C}_{{\theta }_z}+{\omega }_{by}{S}_{{\theta }_z}\\ -{\omega }_{bx}{S}_{{\theta }_z}+{\omega }_{by}{C}_{{\theta }_z}\\ {\dot{\theta }}_z+{\omega }_{bz}\end{array}\right]$$

(1)

where $S_{{\theta }_i}=sin{\theta }_i$, $C_{{\theta }_i}=cos{\theta }_i,i=x,y,z$ and ${({\theta }_i)}_i,i=x,y,z$ represents the rotation matrix of ${\theta }_i$ about the $z_m$-axis, $x_g$-axis and $y_g$-axis, respectively.

Similarly, the angular velocity of the gimbal can be obtained by rotating the motor shaft speed, ${\dot{\theta }}_z$, into the gimbal body frame and adding the rate about the x-axis, so we have:

$${\mathit{\omega }}_{\mathit{g}}=\left[\begin{array}{c}{\omega }_{gx}\\ {\omega }_{gy}\\ {\omega }_{gz}\end{array}\right]=\left[\begin{array}{c}{\dot{\theta }}_x\\ 0\\ 0\end{array}\right]+{({\theta }_x)}_x·{\mathit{\omega }}_{\mathit{m}}=\left[\begin{array}{c}{\dot{\theta }}_x+{\omega }_{bx}{C}_{{\theta }_z}+{\omega }_{by}{S}_{{\theta }_z}\\ -{\omega }_{bx}{C}_{{\theta }_x}{S}_{{\theta }_z}+{\omega }_{by}{C}_{{\theta }_x}{C}_{{\theta }_z}+({\dot{\theta }}_z+{\omega }_{bz})S_{{\theta }_x}\\ {\omega }_{bx}{S}_{{\theta }_x}{S}_{{\theta }_z}-{\omega }_{by}{S}_{{\theta }_x}{C}_{{\theta }_z}+({\dot{\theta }}_z+{\omega }_{bz})C_{{\theta }_x}\end{array}\right]$$

(2)

Finally, rotating the gimbal angular velocities into the rotor body frame and then adding the rate about the y-axis, the angular velocities of the rotor are given as follows:

$$\begin{array}{cc}\hfill {\mathit{\omega }}_{\mathit{r}}& =\left[\begin{array}{c}{\omega }_{rx}\\ {\omega }_{ry}\\ {\omega }_{rz}\end{array}\right]=\left[\begin{array}{c}0\\ {\dot{\theta }}_y\\ 0\end{array}\right]+{\left({\theta }_y\right)}_y·{\mathit{\omega }}_{\mathit{g}}\hfill \\ & =\left[\begin{array}{c}{\dot{\theta }}_x{C}_{{\theta }_y}-{\dot{\theta }}_z{C}_{{\theta }_x}{S}_{{\theta }_y}+\left(C_{{\theta }_y}{C}_{{\theta }_z}-S_{{\theta }_x}{S}_{{\theta }_y}{S}_{{\theta }_z}\right){\omega }_{bx}+\left(C_{{\theta }_y}{S}_{{\theta }_z}+S_{{\theta }_x}{S}_{{\theta }_y}{C}_{{\theta }_z}\right){\omega }_{by}-C_{{\theta }_x}{S}_{{\theta }_y}{\omega }_{bz}\hfill \\ {\dot{\theta }}_z{S}_{{\theta }_x}+{\dot{\theta }}_y-C_{{\theta }_x}{S}_{{\theta }_z}{\omega }_{bx}+C_{{\theta }_x}{C}_{{\theta }_z}{\omega }_{by}+S_{{\theta }_x}{\omega }_{bz}\hfill \\ {\dot{\theta }}_x{S}_{{\theta }_y}+{\dot{\theta }}_z{C}_{{\theta }_x}{C}_{{\theta }_y}+\left(S_{{\theta }_y}{C}_{{\theta }_z}+S_{{\theta }_x}{C}_{{\theta }_y}{S}_{{\theta }_z}\right){\omega }_{bx}+\left(S_{{\theta }_y}{S}_{{\theta }_z}-S_{{\theta }_x}{C}_{{\theta }_y}{C}_{{\theta }_z}\right){\omega }_{by}+C_{{\theta }_x}{C}_{{\theta }_y}{\omega }_{bz}\hfill \end{array}\right]\hfill \end{array}$$

(3)

According to the calculated angular rates of the above different GW bodies, the kinetic energy T of the GW system can be expressed as the generalized rotation speed quadratic forms:

$$T=\frac{1}{2}(\sum _{i=x,y,z}{I}_{mi}{\omega }_{mi}^2+\sum _{i=x,y,z}{I}_{gi}{\omega }_{gi}^2+\sum _{i=x,y,z}{I}_{ri}{\omega }_{ri}^2)$$

(4)

where $I_{ri}$ and $I_{gi}$,$i=x,y,z$ are the moments of inertia of the rotor and gimbal along the corresponding axes, respectively.

Furthermore, the potential energy V is the sum of the potential energy of the inner and outer torsion deformation, which is given by:

$$V=k_x{\theta }_x^2+k_y{\theta }_y^2$$

(5)

Forming the Lagrangian function $L=T-V$ and applying Lagrange’s equations over ${\theta }_x$ and ${\theta }_y$ yield the GW dynamics equations along the transverse axes as follows:

$$\begin{array}{cc}\hfill I_1·{\ddot{\theta }}_x& =-c_x{\dot{\theta }}_x-k_x{\theta }_x-\frac{1}{2}{I}_2{S}_{2{\theta }_x}·{{\dot{\theta }}_z}^2-I_3{S}_{2{\theta }_y}·{\dot{\theta }}_x{\dot{\theta }}_y\hfill \\ & -(I_3{C}_{2{\theta }_y}-I_{ry})C_{{\theta }_x}·{\dot{\theta }}_y{\dot{\theta }}_z+T_{gx}\hfill \\ & +B_1\left(\theta \right){\omega }_{bx}+B_2\left(\theta \right){\omega }_{by}+B_3\left(\theta \right){\dot{\omega }}_{bx}+B_4\left(\theta \right){\dot{\omega }}_{by}\hfill \\ & +B_5\left(\theta \right){\omega }_{bx}^2+B_6\left(\theta \right){\omega }_{by}^2+B_7\left(\theta \right){\omega }_{by}^2\hfill \\ \hfill I_{ry}·{\ddot{\theta }}_y& =-c_y{\dot{\theta }}_y-k_y{\theta }_y-\frac{1}{2}{I}_3{C}_{{\theta }_x}^2{S}_{2{\theta }_y}·{\dot{\theta }}_z^2+\frac{1}{2}{I}_3{S}_{2{\theta }_y}·{\dot{\theta }}_x^2\hfill \\ & +(I_3{C}_{2{\theta }_y}-I_{ry})C_{{\theta }_x}·{\dot{\theta }}_x{\dot{\theta }}_z+T_{gy}\hfill \\ & +D_1\left(\theta \right){\omega }_{bx}+D_2\left(\theta \right){\omega }_{by}+D_3\left(\theta \right){\dot{\omega }}_{bx}+D_4\left(\theta \right){\dot{\omega }}_{by}\hfill \\ & +D_5\left(\theta \right){\omega }_{bx}^2+D_6\left(\theta \right){\omega }_{by}^2+D_7\left(\theta \right){\omega }_{by}^2\hfill \end{array}$$

(6)

where $I_1=I_{gx}+I_{rx}{C}_{{\theta }_y}^2+I_{rz}{S}_{{\theta }_y}^2,I_2=I_{gz}-I_{gy}-I_{ry}+I_{rx}{S}_{{\theta }_y}^2+I_{rz}{C}_{{\theta }_y}^2,I_3=I_{rz}-I_{rx}$ and $B_i$, $D_i$, $i=1,2\cdots ,7$ are nonlinear coefficients in terms of the spacecraft angular rates (${\omega }_{bx},{\omega }_{by}$) as follows:

$$\begin{array}{cc}\hfill B_1=& \left[-I_{re}{S}_{\theta x}{S}_{2{\theta }_y}{C}_{{\theta }_z}+I_1{S}_{{\theta }_z}+I_{2}cos2x_1{S}_{{\theta }_z}\right]{\dot{\theta }}_z\hfill \\ & -\left[I_{re}\left(S_{2{\theta }_y}{C}_{{\theta }_z}+S_{\theta x}{C}_{2{\theta }_y}{S}_{{\theta }_z}\right)-I_{ry}{S}_{\theta x}{S}_{{\theta }_z}\right]{\dot{\theta }}_y\hfill \\ \hfill B_2=& -\left[I_{re}{S}_{\theta x}{S}_{2{\theta }_y}{S}_{{\theta }_z}+I_1{C}_{{\theta }_z}+I_2{C}_{2{\theta }_y}{C}_{{\theta }_z}\right]{\dot{\theta }}_z\hfill \\ & -\left[I_{re}\left(S_{2{\theta }_y}{S}_{{\theta }_z}-S_{\theta x}{C}_{2{\theta }_y}{C}_{{\theta }_z}\right)+I_{ry}{S}_{\theta x}{C}_{{\theta }_z}\right]{\dot{\theta }}_y\hfill \\ \hfill B_3=& -I_1{C}_{{\theta }_z}-\frac{1}{2}{I}_{re}{S}_{\theta x}{S}_{2{\theta }_y}{S}_{{\theta }_z}\hfill \\ \hfill B_4=& -I_1{S}_{{\theta }_z}+\frac{1}{2}{I}_{re}{S}_{\theta x}{S}_{2{\theta }_y}{C}_{{\theta }_z}\hfill \\ \hfill B_5=& \frac{1}{2}{I}_2{S}_{{\theta }_z}^2{S}_{2{\theta }_x}+\frac{1}{4}{I}_{re}{C}_{\theta x}{S}_{2{\theta }_y}{S}_{2{\theta }_z}\hfill \\ \hfill B_6=& \frac{1}{2}{I}_2{C}_{{\theta }_z}^2{S}_{2{\theta }_x}-\frac{1}{4}{I}_{re}{C}_{\theta x}{S}_{2{\theta }_y}{S}_{2{\theta }_z}\hfill \\ \hfill B_7=& -\frac{1}{2}{I}_2{S}_{2{\theta }_z}{S}_{2{\theta }_x}-\frac{1}{2}{I}_{re}{C}_{\theta x}{S}_{2{\theta }_y}{C}_{2{\theta }_z}\hfill \end{array}$$

$$\begin{array}{cc}& D_1=\left[I_{re}\left(S_{2{\theta }_y}{C}_{{\theta }_z}+S_{{\theta }_x}{C}_{2{\theta }_y}{S}_{{\theta }_z}\right)-I_{ry}{S}_{{\theta }_x}{S}_{{\theta }_z}\right]{\dot{\theta }}_x\hfill \\ & +\left[I_{re}\left(C_{{\theta }_x}{C}_{2{\theta }_y}{C}_{{\theta }_z}-\frac{1}{2}{S}_{2{\theta }_x}{S}_{2{\theta }_y}{S}_{{\theta }_z}\right)+I_{ry}{C}_{{\theta }_x}{C}_{{\theta }_z}\right]{\dot{\theta }}_z\hfill \\ & D_2=\left[I_{re}\left(S_{2{\theta }_y}{S}_{{\theta }_z}-S_{{\theta }_x}{C}_{2{\theta }_y}{C}_{{\theta }_z}\right)+I_{ry}{S}_{{\theta }_x}{C}_{{\theta }_z}\right]{\dot{\theta }}_x\hfill \\ & +\left[I_{re}\left(\frac{1}{2}{S}_{2{\theta }_x}{S}_{2{\theta }_y}{C}_{{\theta }_z}+C_{{\theta }_x}{C}_{2{\theta }_y}{S}_{{\theta }_z}\right)+I_{ry}{C}_{{\theta }_x}{S}_{{\theta }_z}\right]{\dot{\theta }}_z\hfill \\ & D_3=I_{ry}{C}_{{\theta }_x}{S}_{{\theta }_z}{D}_4=-I_{ry}{C}_{{\theta }_x}{C}_{{\theta }_z}\hfill \\ & D_5=\frac{1}{2}{I}_{re}\left[S_{{\theta }_x}{C}_{2{\theta }_y}{S}_{2{\theta }_z}+S_{2{\theta }_y}\left(C_{{\theta }_z}^2-S_{{\theta }_x}^2{S}_{{\theta }_z}^2\right)\right]\hfill \\ & D_6=-\frac{1}{2}{I}_{re}\left[S_{{\theta }_x}{C}_{2{\theta }_y}{S}_{2{\theta }_z}-S_{2{\theta }_y}\left(S_{{\theta }_z}^2-S_{{\theta }_x}^2{C}_{{\theta }_z}^2\right)\right]\hfill \\ & D_7=\frac{1}{2}{I}_{re}\left(S_{{\theta }_x}^2{S}_{2{\theta }_y}{S}_{2{\theta }_z}-2S_{{\theta }_x}{C}_{2{\theta }_y}{C}_{2{\theta }_z}+S_{2{\theta }_y}{S}_{2{\theta }_z}\right)\hfill \end{array}$$

Actually, the generalized coordinates $({\theta }_x,{\theta }_y)$ in Equation (6) represent the rotation angles of the inner and outer torsion springs, which cannot be indirectly measured. Therefore, another set of coordinates $({\varphi }_x,{\varphi }_y)$ named “case coordinates” should be defined in the case frame ${\mathcal{F}}_0$, and the case coordinates $({\varphi }_x,{\varphi }_y)$ physically represent the tilt angles of the rotor along the $O{x}_c$-axis and $O{y}_c$-axis, which can be measured by the tilt sensors of GW. Then, the relationships between the generalized coordinates and the case coordinates should be established. The detailed derivation will not be shown in this paper, and the relationships given in the development of DTG are as follows:

$$\begin{array}{cc}\hfill {\theta }_x& ={\varphi }_x{C}_{{\theta }_z}+{\varphi }_y{S}_{{\theta }_z}\hfill \\ \hfill {\theta }_y& =-{\varphi }_x{S}_{{\theta }_z}+{\varphi }_y{C}_{{\theta }_z}\hfill \end{array}$$

(7)

The variables (${\dot{\theta }}_x,{\dot{\theta }}_y$) and (${\ddot{\theta }}_x,{\ddot{\theta }}_y$) in Equation (6) can be calculated by taking the first and second derivatives of Equation (7), then we have:

$$\begin{array}{cc}\hfill {\dot{\theta }}_x& ={\dot{\varphi }}_x{C}_{{\theta }_z}-{\varphi }_x{\dot{\theta }}_z{S}_{{\theta }_z}+{\dot{\varphi }}_y{S}_{{\theta }_z}+{\varphi }_y{\dot{\theta }}_z{C}_{{\theta }_z}\hfill \\ \hfill {\dot{\theta }}_y& =-{\dot{\varphi }}_x{S}_{{\theta }_z}-{\varphi }_x{\dot{\theta }}_z{C}_{{\theta }_z}+{\dot{\varphi }}_y{C}_{{\theta }_z}-{\varphi }_y{\dot{\theta }}_z{S}_{{\theta }_z}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\ddot{\theta }}_x=& {\ddot{\varphi }}_x{C}_{{\theta }_z}-{\dot{\varphi }}_x{\dot{\theta }}_z{S}_{{\theta }_z}-({\dot{\varphi }}_x{\dot{\theta }}_z+{\varphi }_x{\ddot{\theta }}_z)S_{{\theta }_z}-{\varphi }_x{\dot{\theta }}_z^2{C}_{{\theta }_z}\hfill \\ & +{\ddot{\varphi }}_y{S}_{{\theta }_z}+{\dot{\varphi }}_y{\dot{\theta }}_z{C}_{{\theta }_z}+({\dot{\varphi }}_y{\dot{\theta }}_z+{\varphi }_y{\ddot{\theta }}_z)C_{{\theta }_z}-{\varphi }_y{\dot{\theta }}_z^2{S}_{{\theta }_z}\hfill \\ \hfill {\ddot{\theta }}_y=& -{\ddot{\varphi }}_x{S}_{{\theta }_z}-{\dot{\varphi }}_x{\dot{\theta }}_z{C}_{{\theta }_z}-({\dot{\varphi }}_x{\dot{\theta }}_z+{\varphi }_x{\ddot{\theta }}_z)C_{{\theta }_z}+{\varphi }_x{\dot{\theta }}_z^2{S}_{{\theta }_z}\hfill \\ & +{\ddot{\varphi }}_y{C}_{{\theta }_z}-{\dot{\varphi }}_y{\dot{\theta }}_z{S}_{{\theta }_z}-({\dot{\varphi }}_y{\dot{\theta }}_z+{\varphi }_y{\ddot{\theta }}_z)S_{{\theta }_z}-{\varphi }_y{\dot{\theta }}_z^2{C}_{{\theta }_z}\hfill \end{array}$$

(9)

At this point, assume the rotor transverse inertias are equal to each other and set as the value $I_{rt}$. Similarly, the gimbal transverse inertias are also equal to each other and set as the values $I_{gt}$, that is,

$$\begin{array}{cc}\hfill I_{rx}& =I_{ry}=I_{rt}\hfill \\ \hfill I_{gx}& =I_{gy}=I_{gt}\hfill \end{array}$$

(10)

Additionally, the inertia of the moments of the rotor and gimbal along the spin axis are re-assigned the variable names $I_{rs}$ and $I_{gs}$ to more clearly distinguish between spin (sub-subscript “s”) and transverse (sub-subscript “t”).

Taking Equations (7)–(9) into Equation (6) and rearranging the results, the motion of the equations of GW expressed by the case coordinates (${\varphi }_x,{\varphi }_y$) are yielded as follows:

$${\mathit{M}}_c\left(x\right)\ddot{\mathit{x}}+{\mathit{C}}_c\left(\mathit{x}\right)\dot{\mathit{x}}={\mathit{Q}}_c\left(\mathit{x}\right){\mathit{T}}_c+{\mathit{F}}_c\left(\mathit{x},\dot{\mathit{x}}\right)+{\mathit{F}}_{\omega }\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b\right)$$

(11)

where:

$$\begin{array}{cc}& \mathit{x}={\left[\begin{array}{cc}{\varphi }_x& {\varphi }_y\end{array}\right]}^T{\mathit{T}}_c={\left[\begin{array}{cc}{T}_{cx}& T_{cy}\end{array}\right]}^T\hfill \\ & {\mathit{M}}_c\left(\mathit{x}\right)=\left[\begin{array}{cc}{I}_1{C}_{{\theta }_z}& I_1{S}_{{\theta }_z}\\ -I_{rt}{S}_{{\theta }_z}& I_{rt}{C}_{{\theta }_z}\end{array}\right]{\mathit{C}}_c\left(\mathit{x}\right)=\left[\begin{array}{cc}{c}_x{C}_{{\theta }_z}& c_x{S}_{{\theta }_z}\\ -c_y{S}_{{\theta }_z}& c_y{C}_{{\theta }_z}\end{array}\right]\hfill \\ & {\mathit{Q}}_c\left(\mathit{x}\right)=\left[\begin{array}{cc}{C}_{{\theta }_z}& S_{{\theta }_z}\\ -S_{{\theta }_z}{C}_{{\theta }_x}& C_{{\theta }_z}{C}_{{\theta }_x}\end{array}\right]{\mathit{F}}_c\left(\mathit{x},\dot{\mathit{x}}\right)=\left[\begin{array}{c}{f}_{c1}\left(\mathit{x},\dot{\mathit{x}}\right)\\ f_{c2}\left(\mathit{x},\dot{\mathit{x}}\right)\end{array}\right]\hfill \\ & {\mathit{F}}_{\omega }\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b\right)=\left[\begin{array}{c}{f}_{\omega 1}\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b\right)\\ f_{\omega 2}\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b\right)\end{array}\right]\hfill \end{array}$$

$$\begin{array}{cc}& f_{c1}\left(\mathit{x},\dot{\mathit{x}}\right)=-M_{c1}·{\ddot{\theta }}_z-K_{c11}·{\varphi }_x-K_{c12}·{\varphi }_y-C_{c11}·{\dot{\varphi }}_x-C_{c12}·{\dot{\varphi }}_y-\frac{1}{2}{I}_2{S}_{2{\theta }_x}{\dot{\theta }}_z^2-I_{re}{S}_{2{\theta }_y}{\dot{\theta }}_x{\dot{\theta }}_y\hfill \\ & f_{c2}\left(\mathit{x},\dot{\mathit{x}}\right)=-M_{c2}·{\ddot{\theta }}_z-K_{c21}·{\varphi }_x-K_{c22}·{\varphi }_y-C_{c21}·{\dot{\varphi }}_x-C_{c22}·{\dot{\varphi }}_y-\frac{1}{2}{I}_{re}{C}_{{\theta }_x}^2{S}_{2{\theta }_y}{\dot{\theta }}_z^2+\frac{1}{2}{I}_{re}{S}_{2{\theta }_y}{\dot{\theta }}_x^2\hfill \\ & \begin{array}{cc}\hfill f_{\omega 1}\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b\right)& =M_1\left(\mathit{x},\dot{\mathit{x}}\right)·{\dot{\omega }}_{bx}+M_2\left(\mathit{x},\dot{\mathit{x}}\right)·{\dot{\omega }}_{by}+M_3\left(\mathit{x},\dot{\mathit{x}}\right)·{\omega }_{bx}+M_4\left(\mathit{x},\dot{\mathit{x}}\right)·{\omega }_{by}\hfill \end{array}\hfill \\ & \begin{array}{cc}\hfill f_{\omega 2}\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b\right)& =N_1\left(\mathit{x},\dot{\mathit{x}}\right)·{\dot{\omega }}_{bx}+N_2\left(\mathit{x},\dot{\mathit{x}}\right)·{\dot{\omega }}_{by}+N_3\left(\mathit{x},\dot{\mathit{x}}\right)·{\omega }_{bx}+N_4\left(\mathit{x},\dot{\mathit{x}}\right)·{\omega }_{by}\hfill \end{array}\hfill \end{array}$$

where:

$$\begin{array}{cc}& M_{c1}=I_1·\left({\varphi }_y{C}_{{\theta }_z}-{\varphi }_x{S}_{{\theta }_z}\right)+\frac{1}{2}{I}_{re}{C}_{{\theta }_x}{S}_{2{\theta }_y}\hfill \\ & M_{c2}=-I_{rt}\left({\varphi }_x{C}_{{\theta }_z}+{\varphi }_y{S}_{{\theta }_z}-S_{{\theta }_x}\right)\hfill \\ & K_{c11}=-I_1{\dot{\theta }}_z^2{C}_{{\theta }_z}-c_x{\dot{\theta }}_z{S}_{{\theta }_z}+k_x{C}_{{\theta }_z}-I_4{\dot{\theta }}_z{C}_{{\theta }_z}\hfill \\ & K_{c12}=-I_1{\dot{\theta }}_z^2{S}_{{\theta }_z}+c_x{\dot{\theta }}_z{C}_{{\theta }_z}+k_x{S}_{{\theta }_z}-I_4{\dot{\theta }}_z{S}_{{\theta }_z}\hfill \\ & K_{c21}=I_{rt}{\dot{\theta }}_z^2{S}_{{\theta }_z}-c_y{\dot{\theta }}_z{C}_{{\theta }_z}-k_y{S}_{{\theta }_z}+I_4{\dot{\theta }}_z{S}_{{\theta }_z}\hfill \\ & K_{c22}=-I_{rt}{\dot{\theta }}_z^2{C}_{{\theta }_z}-c_y{\dot{\theta }}_z{S}_{{\theta }_z}+k_y{C}_{{\theta }_z}-I_4{\dot{\theta }}_z{C}_{{\theta }_z}\hfill \\ & C_{c11}=-2I_1{\dot{\theta }}_z{S}_{{\theta }_z}-I_4{S}_{{\theta }_z}{C}_{c12}=2I_1{\dot{\theta }}_z{C}_{{\theta }_z}+I_4{C}_{{\theta }_z}\hfill \\ & C_{c21}=-2I_{rt}{\dot{\theta }}_z{C}_{{\theta }_z}-I_4{C}_{{\theta }_z}{C}_{c22}=-2I_{rt}{\dot{\theta }}_z{S}_{{\theta }_z}-I_4{S}_{{\theta }_z}\hfill \\ & I_4=\left(I_{re}cos\left(2{\theta }_y\right)-I_{rt}\right){\dot{\theta }}_z\cos {\theta }_x\hfill \end{array}$$

$M_i\left(\mathit{x},\dot{\mathit{x}}\right),N_i\left(\mathit{x},\dot{\mathit{x}}\right),i=1,\cdots ,4$ are equal to the corresponding expressions $B_i,D_i,i=1,\cdots ,4$; however, the variables ${\theta }_j,{\dot{\theta }}_j,j=x,y$ in $B_i,D_i$ are substituted with Equations (7)–(9).

3. Modeling of Angular Rate Sensing with a Gyrowheel

3.1. Measurement Model at Arbitrary Operating Position Based on Real-Time Lyapunov Linearization

Suppose that a physical system can be generally expressed by the following nonlinear autonomous system equation:

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\right)+\mathit{g}\left(\mathit{x}\right)\mathit{u}$$

(12)

where $\mathit{x}={\left[x_1,x_2,\cdots ,x_n\right]}^T\in {\mathcal{R}}^{n\times 1}$ is the state vector, $\mathit{u}={\left[u_1,u_2,\cdots ,u_m\right]}^T\in {\mathcal{R}}^{m\times 1}$ is the control input vector and $f(\mathit{x})\in {\mathcal{R}}^{n\times 1}$ and $g(\mathit{x})\in {\mathcal{R}}^{n\times 1}$ are vector functions of states.

Assume ${\mathit{x}}_d(t)\in {\mathcal{R}}^{n\times 1}$ is a given reference trajectory whose corresponding reference input is ${\mathit{u}}_d$, then we have:

$${\dot{\mathit{x}}}_d=\mathit{f}({\mathit{x}}_d)+\mathit{g}({\mathit{x}}_d){\mathit{u}}_d$$

(13)

Taking Lyapunov’s linearization [24,25] around the operating points $({\mathit{x}}_d,{\mathit{u}}_d)$, then it yields:

$$\dot{\mathit{x}}={\dot{\mathit{x}}}_d+\mathit{A}\left({\mathit{x}}_d\right)\left(\mathit{x}-{\mathit{x}}_d\right)+\mathit{B}\left({\mathit{x}}_d\right)\left(\mathit{u}-{\mathit{u}}_d\right)$$

(14)

where $\mathit{A}\left({\mathit{x}}_d\right)={\left.\frac{d\mathit{f}}{d\mathit{x}}\right|}_{\mathit{x}={\mathit{x}}_d},\mathit{B}\left({\mathit{x}}_d\right)=\mathit{g}\left({\mathit{x}}_d\right)$.

For the GW system, let ${\mathit{x}}_d,{\dot{\mathit{x}}}_d,{\ddot{\mathit{x}}}_d,{\mathit{\omega }}_{bd},{\dot{\mathit{\omega }}}_{bd}$ be the operating points, which are given by:

$${\mathit{x}}_d^T=\left[\begin{array}{c}{\varphi }_{xd}\\ {\varphi }_{yd}\end{array}\right],{\dot{\mathit{x}}}_d^T=\left[\begin{array}{c}0\\ 0\end{array}\right],{\ddot{\mathit{x}}}_d^T=\left[\begin{array}{c}0\\ 0\end{array}\right],{\mathit{\omega }}_{bd}^T=\left[\begin{array}{c}0\\ 0\end{array}\right],{\dot{\mathit{\omega }}}_{bd}^T=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

(15)

According to Lyapunov’s linearization theory expressed by Equations (13) and (14), the linearized dynamics equations of the GW are given by:

$$\begin{array}{cc}& \frac{\partial \left({\mathit{M}}_c\left(\mathit{x}\right)\ddot{\mathit{x}}\right)}{\partial \ddot{\mathit{x}}}{|}_{\mathit{x}={\mathit{x}}_d}\left(\ddot{\mathit{x}}-{\ddot{\mathit{x}}}_d\right)+\frac{\partial \left({\mathit{M}}_c\left(\mathit{x}\right)\ddot{\mathit{x}}\right)}{\partial \mathit{x}}{|}_{\mathit{x}={\mathit{x}}_d}\left(\mathit{x}-{\mathit{x}}_d\right)\hfill \\ & +\frac{\partial \left({\mathit{C}}_c\left(\mathit{x}\right)\dot{\mathit{x}}\right)}{\partial \dot{\mathit{x}}}{|}_{\mathit{x}={\mathit{x}}_d}\left(\dot{\mathit{x}}-{\dot{\mathit{x}}}_d\right)=\left(\frac{\partial {\mathit{F}}_c}{\partial \mathit{x}}\right){|}_{\begin{array}{c}\hfill \mathit{x}={\mathit{x}}_d\\ \hfill \dot{\mathit{x}}={\dot{\mathit{x}}}_d\end{array}}\left(\mathit{x}-{\mathit{x}}_d\right)\hfill \\ & +\left(\frac{\partial {\mathit{F}}_c}{\partial \dot{\mathit{x}}}\right){|}_{\begin{array}{c}\hfill \mathit{x}={\mathit{x}}_d\\ \hfill \dot{\mathit{x}}={\dot{\mathit{x}}}_d\end{array}}\left(\dot{\mathit{x}}-{\dot{\mathit{x}}}_d\right)+\left(\frac{\partial {\mathit{F}}_{\omega }}{\partial {\mathit{\omega }}_b}\right){|}_{\begin{array}{c}\hfill \mathit{x}={\mathit{x}}_d\\ \hfill \dot{\mathit{x}}={\dot{\mathit{x}}}_d\end{array}}\left({\mathit{\omega }}_b-{\mathit{\omega }}_{bd}\right)\hfill \\ & +\left(\frac{\partial {\mathit{F}}_{\omega }}{\partial {\dot{\mathit{\omega }}}_b}\right){|}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\left({\dot{\mathit{\omega }}}_b-{\dot{\mathit{\omega }}}_{bd}\right)\hfill \\ & +{\mathit{Q}}_c({\mathit{x}}_d)·\left({\mathit{T}}_c-{\mathit{T}}_{cd}\right)+{\mathit{F}}_{h.o.t}\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_{bd},{\dot{\mathit{\omega }}}_{bd}\right)\hfill \end{array}$$

(16)

where:

$$\begin{array}{c}\hfill \frac{\partial \left({\mathit{M}}_c\left(\mathit{x}\right)\ddot{\mathit{x}}\right)}{\partial \ddot{\mathit{x}}}{|}_{\mathit{x}={\mathit{x}}_d}={\left[\begin{array}{cc}{I}_1{C}_{{\theta }_z}& I_1{S}_{{\theta }_z}\\ -I_{rt}{S}_{{\theta }_z}& I_{rt}{C}_{{\theta }_z}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d}{\mathit{Q}}_c\left({\mathit{x}}_d\right)={\left[\begin{array}{cc}{C}_{{\theta }_z}& S_{{\theta }_z}\\ -S_{{\theta }_z}{C}_{{\theta }_x}& C_{{\theta }_z}{C}_{{\theta }_x}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d}\end{array}$$

$$\begin{array}{c}\hfill \frac{\partial \left({\mathit{M}}_c\left(\mathit{x}\right)\ddot{\mathit{x}}\right)}{\partial \mathit{x}}{|}_{\mathit{x}={\mathit{x}}_d}={\left[\begin{array}{cc}-I_{re}{S}_{2{\theta }_y}{S}_{{\theta }_z}\left(C_{{\theta }_z}{\ddot{\varphi }}_x+S_{{\theta }_z}{\ddot{\varphi }}_y\right)& I_{re}{S}_{2{\theta }_y}{C}_{{\theta }_z}\left(C_{{\theta }_z}{\ddot{\varphi }}_x+S_{{\theta }_z}{\ddot{\varphi }}_y\right)\\ 0& 0\end{array}\right]}_{\mathit{x}={\mathit{x}}_d}\end{array}$$

$$\begin{array}{cc}& \left(\frac{\partial {\mathit{F}}_c}{\partial \mathit{x}}\right){|}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}={\left[\begin{array}{cc}\frac{\partial f_{c1}}{\partial {\varphi }_x}& \frac{\partial f_{c1}}{\partial {\varphi }_y}\\ \frac{\partial f_{c2}}{\partial {\varphi }_x}& \frac{\partial f_{c2}}{\partial {\varphi }_y}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\left(\frac{\partial {\mathit{F}}_c}{\partial \dot{\mathit{x}}}\right){|}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}={\left[\begin{array}{cc}\frac{\partial f_{c1}}{\partial {\dot{\varphi }}_x}& \frac{\partial f_{c1}}{\partial {\dot{\varphi }}_y}\\ \frac{\partial f_{c2}}{\partial {\dot{\varphi }}_x}& \frac{\partial f_{c2}}{\partial {\dot{\varphi }}_y}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\hfill \\ & \left(\frac{\partial {\mathit{F}}_{\omega }}{\partial {\mathit{\omega }}_b}\right){|}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}={\left[\begin{array}{cc}\frac{\partial f_{c1}}{\partial {\omega }_{bx}}& \frac{\partial f_{c1}}{\partial {\omega }_{by}}\\ \frac{\partial f_{c2}}{\partial {\omega }_{bx}}& \frac{\partial f_{c2}}{\partial {\omega }_{by}}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\left(\frac{\partial {\mathit{F}}_{\omega }}{\partial {\dot{\mathit{\omega }}}_b}\right){|}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}={\left[\begin{array}{cc}\frac{\partial f_{c1}}{\partial {\dot{\omega }}_{bx}}& \frac{\partial f_{c1}}{\partial {\dot{\omega }}_{by}}\\ \frac{\partial f_{c2}}{\partial {\dot{\omega }}_{bx}}& \frac{\partial f_{c2}}{\partial {\dot{\omega }}_{by}}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\hfill \end{array}$$

${\mathit{F}}_{h.o.t}\left(\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_{bd},{\dot{\mathit{\omega }}}_{bd}\right)$ are the high order terms and will be ignored in the following. For brevity, the elements of the Jacobian matrices of the vectors ${\mathit{F}}_c,{\mathit{F}}_{\omega }$ over the vectors $\mathit{x},\dot{\mathit{x}},{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b$, such as $\frac{\partial f_{ci}}{\partial {\varphi }_j},\frac{\partial f_{ci}}{\partial {\dot{\varphi }}_j},\frac{\partial f_{ci}}{\partial {\omega }_{bj}},\frac{\partial f_{ci}}{\partial {\dot{\omega }}_{bj}},i=1,2,j=x,y$, are also ignored here.

For the equilibrium points of the control input torque vector ${\mathit{T}}_{cd}$, the following constraint condition holds:

$$\begin{array}{cc}\hfill {\mathit{M}}_c({\mathit{x}}_d){\ddot{\mathit{x}}}_d& +{\mathit{C}}_c({\mathit{x}}_d){\dot{\mathit{x}}}_d={\mathit{Q}}_c({\mathit{x}}_d){\mathit{T}}_{cd}+{\mathit{F}}_c({\mathit{x}}_d,{\dot{\mathit{x}}}_d)+{\mathit{F}}_{\omega }({\mathit{x}}_d,{\dot{\mathit{x}}}_d,{\mathit{\omega }}_b,{\dot{\mathit{\omega }}}_b)\hfill \end{array}$$

(17)

Substituting Equation (15) into Equation (17), we have the constraint relationship between the operating state ${\mathit{x}}_d$ and the nominal control input torque ${\mathit{T}}_{cd}$ as follows:

$${\mathit{T}}_{cd}=-{\mathit{Q}}_c^{-1}\left({\mathit{x}}_d\right)·\left({\mathit{F}}_c({\mathit{x}}_d,\mathbf{0})+{\mathit{F}}_{\omega }({\mathit{x}}_d,\mathbf{0},\mathbf{0},\mathbf{0})\right)$$

(18)

Finally, combining Equation (16) with Equation (15) and rearranging the results, the yielded GW Lyapunov linearized equations are given by:

$$\begin{array}{cc}\hfill {\mathit{B}}_c^{\prime }\left({\mathit{x}}_d\right)\left[\begin{array}{c}{\dot{\omega }}_{bx}\\ {\dot{\omega }}_{by}\end{array}\right]& +{\mathit{D}}_c^{\prime }\left({\mathit{x}}_d\right)\left[\begin{array}{c}{\omega }_{bx}\\ {\omega }_{by}\end{array}\right]={\mathit{M}}_c^{\prime }\left({\mathit{x}}_d\right)\left[\begin{array}{c}{\ddot{\varphi }}_x\\ {\ddot{\varphi }}_y\end{array}\right]+{\mathit{C}}_c^{\prime }\left({\mathit{x}}_d\right)\left[\begin{array}{c}{\dot{\varphi }}_x\\ {\dot{\varphi }}_y\end{array}\right]\hfill \\ & +{\mathit{K}}_c^{\prime }\left({\mathit{x}}_d\right)\left[\begin{array}{c}{\varphi }_x-{\varphi }_{xd}\\ {\varphi }_y-{\varphi }_{yd}\end{array}\right]-{\mathit{Q}}_c\left({\mathit{x}}_d\right)\left[\begin{array}{c}{T}_{cx}-T_{cxd}\\ T_{cy}-T_{cyd}\end{array}\right]\hfill \end{array}$$

(19)

where:

$$\begin{array}{cc}\hfill {\mathit{M}}_c^{\prime }\left({\mathit{x}}_d\right)& ={\left[\begin{array}{cc}{I}_1{C}_{{\theta }_z}+{\chi }_1({\mathit{x}}_d)& I_1{S}_{{\theta }_z}+{\chi }_2({\mathit{x}}_d)\\ -I_{rt}{S}_{{\theta }_z}& I_{rt}{C}_{{\theta }_z}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d}\hfill \\ \hfill {\mathit{C}}_c^{\prime }\left({\mathit{x}}_d\right)& ={\left[\begin{array}{cc}{c}_x{C}_{{\theta }_z}-\frac{\partial f_{c1}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\dot{\varphi }}_x}& c_x{S}_{{\theta }_z}-\frac{\partial f_{c1}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\dot{\varphi }}_y}\\ -c_y{S}_{{\theta }_z}-\frac{\partial f_{c2}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\dot{\varphi }}_x}& c_y{C}_{{\theta }_z}-\frac{\partial f_{c2}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\dot{\varphi }}_y}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\hfill \\ \hfill {\mathit{K}}_c^{\prime }\left({\mathit{x}}_d\right)& ={\left[\begin{array}{cc}\frac{\partial f_{c1}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\varphi }_x}& \frac{\partial f_{c1}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\varphi }_y}\\ \frac{\partial f_{c2}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\varphi }_x}& \frac{\partial f_{c2}\left(\mathit{x},\dot{\mathit{x}}\right)}{\partial {\varphi }_y}\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\hfill \\ \hfill {\mathit{B}}_c^{\prime }\left({\mathit{x}}_d\right)& ={\left[\begin{array}{cc}{M}_1\left(\mathit{x},\dot{\mathit{x}}\right)& M_2\left(\mathit{x},\dot{\mathit{x}}\right)\\ N_1\left(\mathit{x},\dot{\mathit{x}}\right)& N_2\left(\mathit{x},\dot{\mathit{x}}\right)\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}{\mathit{D}}_c^{\prime }\left({\mathit{x}}_d\right)={\left[\begin{array}{cc}{M}_3\left(\mathit{x},\dot{\mathit{x}}\right)& M_4\left(\mathit{x},\dot{\mathit{x}}\right)\\ N_3\left(\mathit{x},\dot{\mathit{x}}\right)& N_4\left(\mathit{x},\dot{\mathit{x}}\right)\end{array}\right]}_{\mathit{x}={\mathit{x}}_d,\dot{\mathit{x}}={\dot{\mathit{x}}}_d}\hfill \\ \hfill {\chi }_1({\mathit{x}}_d)& =I_{re}{S}_{2{\theta }_y}·\left(-\frac{1}{2}{S}_{2{\theta }_z}({\varphi }_x-{\varphi }_{xd})-C_{{\theta }_z}^2({\varphi }_y-{\varphi }_{yd})\right)\hfill \\ \hfill {\chi }_2({\mathit{x}}_d)& =I_{re}{S}_{2{\theta }_y}·\left(-S_{2{\theta }_z}^2({\varphi }_x-{\varphi }_{xd})+\frac{1}{2}{S}_{2{\theta }_z}^2({\varphi }_y-{\varphi }_{yd})\right)\hfill \end{array}$$

When the operating points ${\mathit{x}}_d$ are given by the tilt control commands, which are supplied by ACS, the real-time Lyapunov’s linearization of the GW measurement equations can be realized in theory.

3.2. Measurement Model within Small Tilt Angles Based on Complex Quantity Transform

Specially, when the operating points ${\mathit{x}}_d$ are further considered to be set as zero tilt angles, that is, ${\mathit{x}}_d^T=\left[\begin{array}{c}{\varphi }_{xd}\\ {\varphi }_{yd}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$, we have:

$$\begin{array}{cc}\hfill \left[\begin{array}{cc}{I}_t{C}_{{\theta }_z}& I_t{S}_{{\theta }_z}\\ -I_{rt}{S}_{{\theta }_z}& I_{rt}{C}_{{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\ddot{\varphi }}_x\\ {\ddot{\varphi }}_y\end{array}\right]& +\left[\begin{array}{cc}{\varphi }_{d11}& {\varphi }_{d12}\\ {\varphi }_{d21}& {\varphi }_{d22}\end{array}\right]\left[\begin{array}{c}{\dot{\varphi }}_x\\ {\dot{\varphi }}_y\end{array}\right]-\left[\begin{array}{cc}{\varphi }_{11}& {\varphi }_{12}\\ {\varphi }_{21}& {\varphi }_{22}\end{array}\right]\left[\begin{array}{c}{\varphi }_x\\ {\varphi }_y\end{array}\right]=\left[\begin{array}{cc}{C}_{{\theta }_z}& S_{{\theta }_z}\\ -S_{{\theta }_z}& C_{{\theta }_z}\end{array}\right]\left[\begin{array}{c}{T}_{cx}\\ T_{cy}\end{array}\right]\hfill \\ & +\left[\begin{array}{cc}-I_t{C}_{{\theta }_z}& -I_t{S}_{{\theta }_z}\\ I_{rt}{S}_{{\theta }_z}& -I_{rt}{C}_{{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\dot{\omega }}_{bx}\\ {\dot{\omega }}_{by}\end{array}\right]+\left[\begin{array}{cc}{I}_s{\dot{\theta }}_z{S}_{{\theta }_z}& -I_s{\dot{\theta }}_z{C}_{{\theta }_z}\\ I_{rs}{\dot{\theta }}_z{C}_{{\theta }_z}& I_{rs}{\dot{\theta }}_z{S}_{{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\omega }_{bx}\\ {\omega }_{by}\end{array}\right]\hfill \end{array}$$

(20)

where:

$$\begin{array}{cc}\hfill I_t& =I_{rt}+I_{gt}{I}_s=I_{rs}+I_{gs}\hfill \\ \hfill {\varphi }_{d11}& =c_x{C}_{{\theta }_z}-\left(I_{rs}+2I_{gt}\right){\dot{\theta }}_z{S}_{{\theta }_z}{\varphi }_{d12}=c_x{S}_{{\theta }_z}+\left(I_{rs}+2I_{gt}\right){\dot{\theta }}_z{C}_{{\theta }_z}\hfill \\ \hfill {\varphi }_{d21}& =-c_y{S}_{{\theta }_z}-I_{rs}{\dot{\theta }}_z{C}_{{\theta }_z}{\varphi }_{d22}=c_y{C}_{{\theta }_z}-I_{rs}{\dot{\theta }}_z{S}_{{\theta }_z}\hfill \\ \hfill {\varphi }_{11}& =\left(I_{gt}+I_{rs}\right){\ddot{\theta }}_z{S}_{{\theta }_z}+c_x{\dot{\theta }}_z{S}_{{\theta }_z}-k_x{C}_{{\theta }_z}+\left(2I_{gt}-I_{gs}\right){\dot{\theta }}_z^2{C}_{{\theta }_z}\hfill \\ \hfill {\varphi }_{12}& =-\left(I_{rs}+I_{gt}\right){\ddot{\theta }}_z{C}_{{\theta }_z}-c_x{\dot{\theta }}_z{C}_{{\theta }_z}-k_x{S}_{{\theta }_z}+\left(2I_{gt}-I_{gs}\right){\dot{\theta }}_z^2{S}_{{\theta }_z}\hfill \\ \hfill {\varphi }_{21}& =I_{rt}{\ddot{\theta }}_z{C}_{{\theta }_z}+c_y{\dot{\theta }}_z{C}_{{\theta }_z}+k_y{S}_{{\theta }_z}{\varphi }_{22}=I_{rt}{\ddot{\theta }}_z{S}_{{\theta }_z}+c_y{\dot{\theta }}_z{S}_{{\theta }_z}-k_y{C}_{{\theta }_z}\hfill \end{array}$$

Since the two-axis tilt angles of the rotor $({\varphi }_x,{\varphi }_y)$, the two-axis tilt control torques $(T_{cx},T_{cy})$ and the defined two-axis spacecraft angular rates $({\omega }_{bx},{\omega }_{by})$ are perpendicular to each other, respectively, to represent the equation set Equation (20) in the form of one single complex equation [26], the following complex quantities are defined by:

$$\begin{array}{c}\hfill \varphi ={\varphi }_x+j{\varphi }_y{T}_c=T_{cx}+j{T}_{cy}\omega ={\omega }_{bx}+j{\omega }_{by}\end{array}$$

(21)

Then, the case-fixed tilt angles, control torques and spacecraft angular rates are expressed, in terms of the newly-defined complex variables and their complex conjugates, as follows:

$$\left\{\begin{array}{c}{\varphi }_x=\frac{1}{2}\left(\overline{\varphi }+\varphi \right)\hfill \\ {\varphi }_y=\frac{j}{2}\left(\overline{\varphi }-\varphi \right)\hfill \end{array}\right.\left\{\begin{array}{c}{T}_{cx}=\frac{1}{2}\left({\overline{T}}_c+T_c\right)\hfill \\ T_{cy}=\frac{j}{2}\left({\overline{T}}_c-T_c\right)\hfill \end{array}\right.\left\{\begin{array}{c}{\omega }_{bx}=\frac{1}{2}\left({\overline{\omega }}_b+{\omega }_b\right)\hfill \\ {\omega }_{by}=\frac{j}{2}\left({\overline{\omega }}_b-{\omega }_b\right)\hfill \end{array}\right.$$

(22)

The first row of Equation (20) is added to $j=\sqrt{-1}$ times the second row of Equation (20), then substituting Equation (22) into the result and considering the relationship $\prime \prime e^{j{\theta }_z}=cos{\theta }_z+jsin{\theta }_z^{\prime \prime }$, the following single differential equation with complex coefficients is yielded:

$$\begin{array}{cc}& \left(I_{rt}+\frac{1}{2}{I}_{gt}\right)\ddot{\varphi }{e}^{-j{\theta }_z}+\frac{1}{2}{I}_{gt}\ddot{\overline{\varphi }}{e}^{j{\theta }_z}+\left(c_g-j\left(I_{rs}+I_{gt}\right){\dot{\theta }}_z\right)\dot{\varphi }{e}^{-j{\theta }_z}+j{I}_{gt}{\dot{\theta }}_z\dot{\overline{\varphi }}{e}^{j{\theta }_z}\hfill \\ & +\left[\frac{K_x+K_y}{2}-\left(I_{gt}-\frac{1}{2}{I}_{gs}\right){\dot{\theta }}_z^2-j\left(c_g{\dot{\theta }}_z+\frac{I_{gt}+I_{rt}+I_{rs}}{2}{\ddot{\theta }}_z\right)\right]\varphi e^{-j{\theta }_z}\hfill \\ & +\left[\frac{K_x-K_y}{2}-\left(I_{gt}-\frac{1}{2}{I}_{gs}\right){\dot{\theta }}_z^2-j\frac{I_{rt}-I_{gt}-I_{rs}}{2}{\ddot{\theta }}_z\right]\overline{\varphi }{e}^{j{\theta }_z}\hfill \\ & =-\left(I_{rt}+\frac{1}{2}{I}_{gt}\right){\dot{\omega }}_b{e}^{-j{\theta }_z}-\frac{1}{2}{I}_{gt}{\dot{\overline{\omega }}}_b{e}^{j{\theta }_z}+j\left(I_{rs}+\frac{1}{2}{I}_{gs}\right){\dot{\theta }}_z{\omega }_b{e}^{-j{\theta }_z}-j\frac{1}{2}{I}_{gs}{\dot{\theta }}_z{\overline{\omega }}_b{e}^{j{\theta }_z}+T_c{e}^{-j{\theta }_z}\hfill \end{array}$$

(23)

Rearranging Equation (23), we have:

$$\begin{array}{c}\hfill e^{-j{\theta }_z}\left\{\begin{array}{cc}& \ddot{\varphi }+\left(\frac{2c_g}{2I_{rt}+I_{gt}}-j\left(\frac{2I_{rs}}{2I_{rt}+I_{gt}}+\frac{2I_{gt}}{2I_{rt}+I_{gt}}\right){\dot{\theta }}_z\right)\dot{\varphi }\hfill \\ & +\left[\frac{K_x+K_y-2\left(I_{gt}-\frac{1}{2}{I}_{gs}\right){\dot{\theta }}_z^2}{2I_{rt}+I_{gt}}-j\left(\frac{2c_g}{2I_{rt}+I_{gt}}{\dot{\theta }}_z+\frac{I_{gt}+I_{rt}+I_{rs}}{2I_{rt}+I_{gt}}{\ddot{\theta }}_z\right)\right]\varphi \hfill \\ & +{\dot{\omega }}_b-j\left(\frac{2I_{rs}}{2I_{rt}+I_{gt}}+\frac{I_{gs}}{2I_{rt}+I_{gt}}\right){\dot{\theta }}_z{\omega }_b\hfill \\ & +e^{j2{\theta }_z}\left\{\begin{array}{cc}& \frac{I_{gt}}{2I_{rt}+I_{gt}}\ddot{\overline{\varphi }}+j\frac{2I_{gt}}{2I_{rt}+I_{gt}}{\dot{\theta }}_z\dot{\overline{\varphi }}\hfill \\ & +\left[\frac{K_x-K_y-2\left(I_{gt}-\frac{1}{2}{I}_{gs}\right){\dot{\theta }}_z^2}{2I_{rt}+I_{gt}}-j\frac{I_{rt}-I_{gt}-I_{rs}}{2I_{rt}+I_{gt}}{\ddot{\theta }}_z\right]\overline{\varphi }\hfill \\ & +\frac{I_{gt}}{2I_{rt}+I_{gt}}{\dot{\overline{\omega }}}_b+j\frac{I_{gs}}{2I_{rt}+I_{gt}}{\dot{\theta }}_z{\overline{\omega }}_b\hfill \end{array}\right\}\hfill \end{array}\right\}=e^{-j{\theta }_z}·\frac{T_c}{I_{rt}+\frac{1}{2}{I}_{gt}}\end{array}$$

(24)

For conciseness, the following quantities are defined by:

$$\begin{array}{cc}\hfill c_g^{\prime }& =\frac{c_g}{2I_{rt}+I_{gt}}{\omega }_n=\frac{2I_{rs}}{2I_{rt}+I_{gt}}{\dot{\theta }}_z\gamma =\frac{I_{gt}}{2I_{rt}+I_{gt}}{\gamma }_2=\frac{I_{gs}}{2I_{rt}+I_{gt}}\hfill \\ \hfill \alpha & =\frac{K_x+K_y-2\left(I_{gt}-\frac{1}{2}{I}_{gs}\right){\dot{\theta }}_z^2}{2I_{rt}+I_{gt}}\beta =\frac{K_x-K_y-2\left(I_{gt}-\frac{1}{2}{I}_{gs}\right){\dot{\theta }}_z^2}{2I_{rt}+I_{gt}}\hfill \\ \hfill I_{rt}^{\prime }& =I_{rt}+\frac{1}{2}{I}_{gt}{J}_p=\frac{I_{gt}+I_{rt}+I_{rs}}{2I_{rt}+I_{gt}}{J}_m=\frac{I_{rt}-I_{gt}-I_{rs}}{2I_{rt}+I_{gt}}\hfill \end{array}$$

and taking the newly-defined quantities into Equation (24) yields:

$$\begin{array}{c}\hfill e^{-j{\theta }_z}\left\{\begin{array}{cc}& \ddot{\varphi }+\left(2c_g^{\prime }-j\left({\omega }_n+2\gamma {\dot{\theta }}_z\right)\right)\dot{\varphi }\hfill \\ & +\left[\alpha -j\left(2C_g^{\prime }{\dot{\theta }}_z+J_p{\ddot{\theta }}_z\right)\right]\varphi +{\dot{\omega }}_b-j\left({\omega }_n+{\gamma }_2{\dot{\theta }}_z\right){\omega }_b\hfill \\ & +e^{j2{\theta }_z}\left\{\gamma \ddot{\overline{\varphi }}+j2\gamma {\dot{\theta }}_z\dot{\overline{\varphi }}+\left(\beta -j{J}_m{\ddot{\theta }}_z\right)\overline{\varphi }+\gamma {\dot{\overline{\omega }}}_b+j{\gamma }_2{\dot{\theta }}_z{\overline{\omega }}_b\right\}\hfill \end{array}\right\}=\frac{T_c}{I_{rt}^{\prime }}·e^{-j{\theta }_z}\end{array}$$

(25)

Eliminating the factor $e^{-j{\theta }_z}$ in Equation (25), we obtain the more concise form of the GW linearization equations at zero tilt angles given by:

$$\begin{array}{cc}& \ddot{\varphi }+\left(2c_g^{\prime }-j\left({\omega }_n+2\gamma {\dot{\theta }}_z\right)\right)\dot{\varphi }+\left[\alpha -j\left(2c_g^{\prime }{\dot{\theta }}_z+J_p{\ddot{\theta }}_z\right)\right]\varphi +{\dot{\omega }}_b-j\left({\omega }_n+{\gamma }_2{\dot{\theta }}_z\right){\omega }_b\hfill \\ & +e^{j2{\theta }_z}\left[\gamma \ddot{\overline{\varphi }}+j2\gamma {\dot{\theta }}_z\dot{\overline{\varphi }}+\left(\beta -j{J}_m{\ddot{\theta }}_z\right)\overline{\varphi }+\gamma {\dot{\overline{\omega }}}_b+j{\gamma }_2{\dot{\theta }}_z{\overline{\omega }}_b\right]=\frac{T_c}{I_{rt}^{\prime }}\hfill \end{array}$$

(26)

Finally, substituting Equation (21) into Equation (26), substituting $e^{2{\theta }_z}$ with $C_{2{\theta }_z}+j{S}_{2{\theta }_z}$ and restoring the results expressed by the complex quantity into real-value equations, we obtain the GW linearization equations at zero tilt angles represented in the real-value form:

$$\begin{array}{cc}& \left[\begin{array}{cc}1+\gamma C_{2{\theta }_z}& \gamma S_{2{\theta }_z}\\ \gamma S_{2{\theta }_z}& 1-\gamma C_{2{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\ddot{\varphi }}_x\\ {\ddot{\varphi }}_y\end{array}\right]+\left[\begin{array}{cc}2c_g^{\prime }-2\gamma {\dot{\theta }}_z{S}_{2{\theta }_z}& \left({\omega }_n+2\gamma {\dot{\theta }}_z\right)+2\gamma {\dot{\theta }}_z{C}_{2{\theta }_z}\\ -\left({\omega }_n+2\gamma {\dot{\theta }}_z\right)+2\gamma {\dot{\theta }}_z{C}_{2{\theta }_z}& 2c_g^{\prime }+2\gamma {\dot{\theta }}_z{S}_{2{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\dot{\varphi }}_x\\ {\dot{\varphi }}_y\end{array}\right]\hfill \\ & +\left[\begin{array}{cc}\alpha +\beta C_{2{\theta }_z}& 2c_g^{\prime }{\dot{\theta }}_z+\beta S_{2{\theta }_z}+\left(J_p-J_m{C}_{2{\theta }_z}+J_m{S}_{2{\theta }_z}\right){\ddot{\theta }}_z\\ -2c_g^{\prime }{\dot{\theta }}_z+\beta S_{2{\theta }_z}-\left(J_p+J_m{C}_{2{\theta }_z}-J_m{S}_{2{\theta }_z}\right){\ddot{\theta }}_z& \alpha -\beta C_{2{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\varphi }_x\\ {\varphi }_y\end{array}\right]\hfill \\ & +\left[\begin{array}{cc}1+\gamma C_{2{\theta }_z}& \gamma S_{2{\theta }_z}\\ \gamma S_{2{\theta }_z}& 1-\gamma C_{2{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\dot{\omega }}_{bx}\\ {\dot{\omega }}_{by}\end{array}\right]+\left[\begin{array}{cc}-{\gamma }_2{\dot{\theta }}_z{S}_{2{\theta }_z}& \left({\omega }_n+{\gamma }_2{\dot{\theta }}_z\right)+{\gamma }_2{\dot{\theta }}_z{C}_{2{\theta }_z}\\ -\left({\omega }_n+{\gamma }_2{\dot{\theta }}_z\right)+{\gamma }_2{\dot{\theta }}_z{C}_{2{\theta }_z}& {\gamma }_2{\dot{\theta }}_z{S}_{2{\theta }_z}\end{array}\right]\left[\begin{array}{c}{\omega }_{bx}\\ {\omega }_{by}\end{array}\right]=\frac{1}{I_{rt}^{{}^{\prime }}}\left[\begin{array}{c}{T}_{cx}\\ T_{cy}\end{array}\right]\hfill \end{array}$$

(27)

From the left side of Equation (27), each element of the coefficient matrices contains the twice periodic components about the motor spin speed ${\dot{\theta }}_z$. Since the periodic terms have no effects on the measurement accuracy of the spacecraft angular rates (${\omega }_{bx},{\omega }_{by}$) and ignoring the periodic terms will save the unnecessary sensor for measuring the motor rotation angle ${\theta }_z$, therefore, the periodic terms are ignored with the model simplifying to the following form:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\ddot{\varphi }}_x\\ {\ddot{\varphi }}_y\end{array}\right]+& \left[\begin{array}{cc}2c_g^{\prime }& \left({\omega }_n+2\gamma {\dot{\theta }}_z\right)\\ -\left({\omega }_n+2\gamma {\dot{\theta }}_z\right)& 2c_g^{\prime }\end{array}\right]\left[\begin{array}{c}{\dot{\varphi }}_x\\ {\dot{\varphi }}_y\end{array}\right]+\left[\begin{array}{cc}\alpha & \left(2c_g^{\prime }{\dot{\theta }}_z+J_p{\ddot{\theta }}_z\right)\\ -\left(2c_g^{\prime }{\dot{\theta }}_z+J_p{\ddot{\theta }}_z\right)& \alpha \end{array}\right]\left[\begin{array}{c}{\varphi }_x\\ {\varphi }_y\end{array}\right]\hfill \\ & +\left[\begin{array}{c}{\dot{\omega }}_{bx}\\ {\dot{\omega }}_{by}\end{array}\right]+\left[\begin{array}{cc}0& \left({\omega }_n+{\gamma }_2{\dot{\theta }}_z\right)\\ -\left({\omega }_n+{\gamma }_2{\dot{\theta }}_z\right)& 0\end{array}\right]\left[\begin{array}{c}{\omega }_{bx}\\ {\omega }_{by}\end{array}\right]=\frac{1}{I_{rt}^{\prime }}\left[\begin{array}{c}{T}_{cx}\\ T_{cy}\end{array}\right]\hfill \end{array}$$

(28)

4. Analysis of the Measurement Schemes of Spacecraft Angular Rates with the GW

In order to analyze the proposed measurement equations, the following two definitions are given:

Definition 1. 

If the tilt angle of the GW rotor works at a certain tilt position and remains unchanged as the spacecraft angular rates are measured with the GW, we call this situation “$StaticMeasurement$”.

Definition 2. 

In contrast to Definition 1, when the tilt angular velocity of the GW rotor is nonzero as the spacecraft angular rates are measured with the GW, we call this situation “$DynamicMeasurement$”.

From the derived measurement Equations (20) and (28), three pieces of valuable information are summarized as follows:

(1)

Compared to Equation (20), the forms of linearization Equation (28) at zero tilt angles are obviously more concise. More importantly, after ignoring the twice periodic components in Equation (27), there exist no terms about motor rotation angle ${\theta }_z$ in Equation (28), which means that it is more convenient for the spacecraft angular rate sensing, since the special sensor for measuring the motor rotation angle in real time is no longer needed in this kind of situation.

(2)

Equation (28) is obtained by substituting the operating points at zero tilt angles into Equation (20), which are suitable for the arbitrary operating state of the rotor. Therefore, combining the analysis (1), in the smaller tilt range of the rotor, Equation (28) is utilized to realize the spacecraft angular rate sensing with GW. However, when the tilt angles of the rotor become larger, the measurement accuracy with Equation (28) cannot meet the indicator requirement without any error compensation. The applications of the real-time linearization measurement Equation (20) and small tilt measurement Equation (28) are analyzed in the following figure.

In Figure 3, considering GW in the application scenario of static measurement, when the rotor works at large tilt angles, the measurement errors due to the linearization at zero tilt angles are significantly correlated to the tilt angles of the rotor, which can be modeled and compensated in the measurement Equation (28) by the methods of polynomial fitting, B-spline functions or table lookups based on the calibration data. Since the rotor needs to keep still for realizing the spacecraft angular rate sensing in this situation, the radial control torques cannot be supplied by the GW, which means that with Equation (28) as the measurement equations of the GW, it has two kinds of working modes: the radial torque outputs mode and the spacecraft angular sensing mode, and these two modes cannot be realized at the same time, we can name the two working modes as $\prime \prime time$-$sharingmultiplexin{g}^{\prime \prime }$. As previously mentioned, the greatest advantages of this scheme are that there is no need to set the rotary transformer for measuring the motor rotation angle ${\theta }_z$, and the regular errors like the linearization at zero tilt angles can be more easily calibrated by ground experiments. However, this scheme has an obvious drawback that GW cannot realize the radial torque outputs and angular rate sense simultaneously, so that the application of GW is limited in the area of spacecraft.

Whereas in the state of dynamic measurement, since the measurement errors are not only relevant to the tilt angles of the rotor, but also to the tilt angular velocities and accelerations, it is difficult for Equation (28) to compensate the linearization errors by polynomial fitting or B-spline functions due to the huge workload. Therefore, in the case of dynamic measurement, the real-time Lyapunov’s linearization measurement Equation (20) need to be applied to sense the spacecraft angular rates, which can avoid the errors of the small tilt linearization Equation (28) at larger tilt angles. Most important of all, gyroscopic moments can be generated at the same time because of the existence of the angular momentum and the tilt angular velocity of the rotor in dynamic measurement. However, the prices are that a rotary transformer must be assembled in the GW system to measure the motor rotation angle ${\theta }_z$, and the real-time disturbance estimation methods like [27,28] should be further developed to compensate regular mechanical errors for Equation (20), which will not be discussed in detail in this paper.

(3)

Since the effects of the motor spin acceleration are considered in Equation (28), while the motor spin speed is being changed to control the spacecraft attitude along the spin axis, the angular rate sensing can be more accurately realized by Equation (28) in small tilt ranges.

5. Simulations

5.1. Simulation Platform

To demonstrate the effectiveness of the proposed angular rate measurement approach, a simulation platform as in Figure 4 is built. The simulation platform is divided into two layers: the torque output layer and the measurement layer. The former contains the motor control loop, the X-axis and the Y-axis tilt control loops. The variables ${\dot{\theta }}_{zd},{\varphi }_{xd},{\varphi }_{yd}$ are the command inputs of the GW control loops. The variables ${\dot{\theta }}_z,{\varphi }_x,{\varphi }_y$ are the corresponding measurable outputs of the above three control loops, respectively. The variables $T_{cz},T_{cx},T_{cy}$ are control torques of the above three control loops, where $T_{cz}$ is used to control the motor rotation speed, and the torques $T_{cx},T_{cy}$ make the rotor tilt along radial directions. The variables $T_z,T_x,T_y$ acting on the spacecraft block are the three-axis control torques generated by the GW system. The sensor block includes the X-axis, Y-axis tilt angle sensors, the currents of torquer coils, motor rotation speed and angle sensors, which can directly measure the variables ${\varphi }_x,{\varphi }_y,T_{cx},T_{cy},{\theta }_z,{\dot{\theta }}_z$, respectively.

The measurement layer is designed to realize the angular rate sense of the spacecraft with the proposed measurement methods, which contains the real-time Lyapunov linearization Equation (20) and the small tilt linearization Equation (28). The tilt command inputs $({\varphi }_{xd},{\varphi }_{yd})$ as the desired operating points are utilized to realize the real-time Lyapunov linearization measurement. The outputs of the sensor block are input to both measurement equations. As previously mentioned, the disturbance estimation (refer to [27,28]) should be studied further for Equation (20) in engineering. The regular error compensation is investigated for Equation (28) for compensating the errors due to zero tilt angle linearization.

The spacecraft angular rates ${\omega }_{bx}^{\prime },{\omega }_{by}^{\prime }$ can be obtained from the proposed measurement equations; however, unlike the zero position linearization Equation (28), the twice periodic components are hardly separated and ignored from the real-time Lyapunov linearization Equation (20). Therefore, if the angular rates are sensed by Equation (20), the twice periodic components in the results should be filtered out by the notch filters $G(s)$, which can be designed as follow:

$$G(s)=\frac{s^2+2\pi as+{(2\pi f)}^2}{s^2+2\pi bs+{(2\pi f)}^2}$$

where f is the center frequency of the notch filter and given by $f=2f_{motor}=50$ Hz; here, $f_{motor}$ represents the spin frequency of the motor; the parameters a and b are given by $0.1$ and 60, respectively.

The key parameters in the simulation are given in

Table 1. Gyrowheel and spacecraft design parameter in simulation.

Parameter Name	Value
Rotor transverse-axis inertia $I_{rt}$	$3.458\times {10}^{-3}$ kg · m${}^2$
Rotor spin-axis inertia $I_{rs}$	$6.402\times {10}^{-3}$ kg · m${}^2$
Gimbal transverse-axis inertia $I_{gt}$	$1.2758\times {10}^{-5}$ kg · m${}^2$
Gimbal spin-axis inertia $I_{gs}$	$1.8047\times {10}^{-5}$ kg · m${}^2$
Torsion spring stiffness $k_x,k_y$	$0.092\mathrm{Nm}/\mathrm{rad}$
Torsion spring damping $c_x,c_y$	0 Nm/(rad/s)
Operating tilt range of rotor ${\varphi }_x,{\varphi }_y$	$0°\le {\varphi }_x,{\varphi }_y\le 5°$
Spacecraft inertia ${\mathit{I}}_{\mathit{sat}}$	$diag\left[\begin{array}{ccc}50& 50& 50\end{array}\right]$ kg · m${}^2$

. Besides, the initial attitude angular velocity of the spacecraft ${\mathit{\omega }}_b=[{\omega }_{bx};{\omega }_{by};{\omega }_{bz}]=[0.001;0.001;0]$ rad/s, and ideally, we assume there exists no disturbance torque acting on the spacecraft, but the control torques from GW.

In the following simulation, the static and dynamic measurements are studied, respectively. In both cases, the proposed measurement Equations (20) and (28) will be analyzed.

5.2. Static Measurement Validation

Considering the situation of the static measurement, which means that there is no radial control torque outputs in this case, the proposed measurement Equations (20) and (28) are applied to measure the spacecraft angular rates at the operating tilt angle range. As an example, when the constant motor speed of $157.04$ rad/s is given and the X-axis tilt control commands of the rotor (${\varphi }_{xd}$) are given by $0,0.1,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5$ in turn, the change curves of the measurement errors of Equations (20) and (28) with the tilt angles are shown in Figure 5.

From Figure 5, we find:

• The measurement accuracy decreases with the increase of the tilt angle for both the real-time linearization equation and the zero position linearization equation; however, due to the linearization operating points at zero tilt angles, the measurement accuracy of the latter deteriorates more seriously than the former.
• The absolute measurement accuracy of the real-time linearization Equation (20) within operating ranges meets the adequate accuracy requirement. Therefore, Equation (20) can be applied to the operating state at arbitrary tilt angles, but if there is not any linearization error compensation for Equation (28), the zero position linearization Equation (28) is more suitable for the small tilt angle of the rotor, rather than the larger tilt angle of the rotor (for example, see Figure 6 and Figure 7.

Further, considering the measurement Equation (28) without the need for sensing the motor rotation angle ${\theta }_z$ and that the fault-tolerant performance of the measurement Equation (28) is more excellent than Equation (20), if the measurement errors of Equation (28) can be compensated, it will be a better choice in the time-sharing multiplexing scheme. Since the linearization errors vary regularly with the tilt angle according to the error curve in Figure 5, the error compensation expressions of the measured spacecraft angular rates (${\omega }_{bxe},{\omega }_{bye}$) are established for Equation (28) based on the simulation error data with polynomial fitting functions, which are given by:

$$\begin{array}{cc}\hfill {\omega }_{bxe}=& -1.447\times {10}^{-7}·{\varphi }_{xd}^5+1.379\times {10}^{-6}·{\varphi }_{xd}^4-6.801\times {10}^{-6}·{\varphi }_{xd}^3\hfill \\ & +1.080\times {10}^{-5}·{\varphi }_{xd}^2-8.840\times {10}^{-6}·{\varphi }_{xd}+5.349\times {10}^{-7}\hfill \\ \hfill {\omega }_{bye}=& -4.355\times {10}^{-7}·{\varphi }_{xd}^5+1.501\times {10}^{-6}·{\varphi }_{xd}^4+2.008\times {10}^{-4}·{\varphi }_{xd}^3\hfill \\ & +6.730\times {10}^{-6}·{\varphi }_{xd}^2-8.229\times {10}^{-6}·{\varphi }_{xd}+8.113\times {10}^{-8}\hfill \end{array}$$

(29)

With Equation (29) as the compensation equations of the linearization errors, that is ${\omega }_{bx\_new}={\omega }_{bx}+{\omega }_{bxe},{\omega }_{by\_new}={\omega }_{by}+{\omega }_{bye}$, then the measurement error curves with the compensated small tilt linearization equations are also shown in Figure 5 and expressed by the green dotted line and squire markers. From this figure, the measurement accuracy with the measurement Equation (28) after the polynomial compensation is improved significantly, even superior to the real-time linearization measurement Equation (20) in this case. Therefore, when the sharing time multiplexing measurement scheme is adopted, compared to Equation (20), the measurement Equation (28) with regular error compensation components is an ideal method of spacecraft angular rate sensing.

Especially, the response curves of the X-axis tilt control commands given by $0.1$° and $4.8$° are shown as Figure 6 and Figure 7, respectively. From Figure 6, when the rotor is working at the tilt angle of $0.1$°, the measurement accuracy of the zero position linearization Equations (28) with compensation terms (${\omega }_{bxe},{\omega }_{bye}$) is improved to precede ${10}^{-7}$ rad/s compared to the ones without any compensation term. At this point, the real-time Lyapunov linearization measurement can also achieve the perfect measurement performance of the spacecraft angular rates. When the tilt angle of the rotor is increased to $4.8$° in Figure 7, the measurement errors of the zero position linearization Equation (28) with compensation terms (${\omega }_{bxe},{\omega }_{bye}$) are far less than ${10}^{-6}$ rad/s. The real-time Lyapunov linearization measurement accuracy can arrive at ${10}^{-6}$ rad/s. However, the measurement accuracy of the zero position linearization Equation (28) without any compensation term deteriorates seriously, so that the uncompensated zero position linearization equations are not suitable for the angular rate sensing at the large tilt angle. However, since the motor spin angle sensor needs to be introduced to realize the measurement of the motor spin angle when the real-time Lyapunov linearization measurement equations are used to realize the angular rate sensing, it will make the measurement realization more complex in the time-sharing multiplexing scheme.

Comparing Figure 6 to Figure 7, the foregoing analysis for static measurement can be validated; the most important and useful information is given: zero tilt angle linearization Equation (28) with the regular error compensation components is the best choice for the time-sharing multiplexing scheme of GW.

5.3. Dynamic Measurement Validation

In this subsection, we investigate the dynamic measurement of the spacecraft angular rates with the GW, in which case, the control command of the X-axis tilt angle ${\varphi }_{xd}$ keeps as time varying and is given by ${\varphi }_{xd}=4sin(2\pi ·0.05t)$°, so that the Y-axis attitude angle of the spacecraft can be driven by the gyroscopic moment from the GW. The response curves of the dynamic measurement of the spacecraft angular rates with the GW are shown in Figure 8.

Similarly with the previous static measurement case, in Figure 8, the X-axis and Y-axis angular rates can be sensed accurately by the real-time Lyapunov’s linearization method using GW, and the measurement equations linearized at zero tilt position no longer satisfy the measurement requirement. Moreover, since the linearization errors depend on the tilt angular rate of the rotor, it is difficult to establish the compensation polynomial equations by simulation data in dynamic measurement. Therefore, even if in the dynamic measurement case, which means that the gyroscopic moments are produced to control the spacecraft attitudes by tilting the rotor of GW along the radial direction, the proposed innovative measurement Equation (20) based on real-time Lyapunov’s linearization is also an effective way.

6. Conclusions

The GW can not only realize the function of the two-dimensional angular rate sensing, but also the three-dimensional torque output at the same time. The angular rate sensing with GW could be an effective supplement to the conventional ACS rate gyro configurations.

To achieve the above-mentioned goal, two principal contributions to the angular rate measurement with GW are made in this paper:

(1)

By combining the real-time Lyapunov’s linearization with the complex quantity transform, two different measurement models of the spacecraft angular rate sensing with GW are established; we named them as the “real-time Lyapunov linearization measurement model” and the “small tilt measurement model”, respectively.

(2)

For both established measurement models of GW, two different application schemes are proposed: “time-sharing multiplexing of actuator and sensor” and “simultaneous realization of actuator and sensor”. Additionally, the advantage and disadvantage of the two measurement schemes is also presented.

Limited to the experiment conditions, simulations are performed instead of the validation experiment with the GW prototype in this paper. By combining the ground calibration, the proposed measurement method needs to be investigated further in the future.