The Angular Momentum Unloading of the Asymmetric GEO Satellite by Using Electric Propulsion with a Mechanical Arm

Abstract

A high-precision attitude control satellite uses an angular momentum exchange device such as a flywheel or a control moment gyro as the actuator for attitude stability control. Once the accumulation of angular momentum exceeds the upper limit of the angular momentum exchange device, the satellite will lose its attitude control ability. Therefore, it is necessary to unload the angular momentum exchange device to ensure the attitude control ability of the satellite platform. The angular momentum accumulation of GEO(Geosynchronous Orbit, GEO) satellites with asymmetric structure can reach 40 Nms per day, and the accumulation speed is more than 20 times that of GEO satellites with symmetrical structure. Therefore, it is necessary to carry out angular momentum unloading for GEO satellites with asymmetric structure every day. The previous method of angular momentum unloading using electric propulsion has weak unloading capacity, which is not suitable for angular momentum unloading of asymmetric satellites. This paper presents a method of angular momentum unloading using a four-joint mechanical arm plus an electric thruster. Large angular momentum unloading with near-zero burn-up can be achieved through the thrust generated by station keeping. In addition, the problem of attitude and orbit coupling control can be solved by controlling the thrust direction of the electric thruster with a mechanical arm.

1. Introduction

A high-precision attitude control satellite uses an angular momentum exchange device such as a flywheel or a control moment gyro as the actuator for attitude stability control. Once the angular momentum accumulation of the satellite exceeds the upper limit of the capability of the angular momentum exchange device, the satellite will lose its attitude control ability, so the angular momentum exchange device needs to be unloaded to ensure the attitude control ability of the satellite platform. Asymmetrical GEO satellites are generally equipped with single-wing sails or carry large asymmetric loads. The angular momentum accumulated on a surface of such a satellite can reach 40 Nms per day, which is more than 20 times faster than that of a symmetrical GEO satellite. Therefore, it is necessary to provide a large unloading force every day to achieve the unloading of the angular momentum of such a satellite.

At present, chemical propulsion is commonly used in engineering to unload angular momentum. However, the specific impulse of chemical propellants is low, and more propellants are consumed to unload the corresponding angular momentum. The specific impulse of EP (Electric Propulsion, EP) is 5~10 times that of chemical propulsion, and the propellant consumption of the same speed increment is only about 10~20% of chemical propulsion. Therefore, the use of EP can increase the life of satellites and help meet the needs of long life of high-orbit satellite platforms. In addition, the EP is smaller than the chemical propulsion. The use of EP for orbit and attitude control and angular momentum unloading has less disturbance torque to the satellite, which is helpful to realize the high-precision attitude and orbit control of the satellite. Therefore, GEO satellites with all-EP have become a trend.

In recent years, many scholars have studied the use of EP for angular momentum unloading. Walsh [1] introduced the GEO satellite unloading scheme using four-vector EP of a Boeing 702SP satellite. The angular momentum unloading of GEO satellites is realized by slightly deflecting the EP vector to form the unloading torque. By iterating the deflection angle of the EP vector several times, the unloading accuracy of angular momentum is further improved. Ma [2] et al. adopted the scheme of four-vector EP of a Boeing 702SP platform. The angular deflection algorithm of the EP for angular momentum unloading in this scheme is optimized. Although the accuracy of angular momentum unloading is slightly decreased, the calculation amount of angular momentum unloading is significantly reduced, which is conducive to the realization of autonomous angular momentum unloading on board. However, the four-vector EP scheme is not suitable for large angular momentum unloading. Weiss [3,4] introduced the angular momentum unloading scheme of GEO satellites using six-fixed EP, which has a large angular momentum unloading capacity. Li et al. [5] adopted the scheme of a double four-degrees-of-freedom robotic arm plus an EP. The possibility of angular momentum unloading by using the propulsion generated by station keeping (SK) of GEO satellites is discussed.

The above scheme of angular momentum unloading by EP mainly has the following problems [6,7,8,9,10,11,12,13,14,15]:

(1)

The angular momentum unloading scheme of four-vector EP has the characteristics of small deflection angle of electric thrust, weak position and attitude adjustment ability, small angular momentum unloading ability, etc., which cannot meet the requirements of angular momentum unloading of asymmetric structure satellites;

(2)

Although the angular momentum unloading scheme with fixed EP can achieve large angular momentum unloading, the thrust of angular momentum unloading cannot be fully used for SK of GEO satellites, resulting in fuel waste. The fixed installation of electric thrusters cannot realize the active adjustment of position and attitude, and the number of electric thrusters to be configured is large, which is difficult to be deployed on GEO satellites.

In the literature [16,17,18], we have studied the NSSK (north–south station keeping, NSSK) and EWSK (east–west station keeping, EWSK) of GEO satellites by using a mechanical arm and EP. This paper uses the mechanical arm to adjust the position and attitude of thrust generated during SK, and the three-axis large angular momentum unloading of asymmetrical satellite is realized. The angular momentum unloading algorithm studied in this paper has the following three innovations:

(1)

In the past, 2- to 4-vector electric thrusters were needed to achieve small angular momentum unloading. In this paper, only one EP with a mechanical arm is needed to achieve large angular momentum unloading, which can reduce the hardware configuration on board;

(2)

In the past, angular momentum unloading required ignition in 2 to 4 different orbital arcs. In this paper, the angular momentum unloading can be achieved only by firing in one orbit arc, which reduces the number of ignition times.

(3)

In the process of attitude maneuvering, SK and angular momentum unloading can be carried out at the same time, which can realize complex attitude and orbit coupling control.

2. Materials and Methods

2.1. Establishment of Coordinate System

2.1.1. Satellite Orbital Coordinate System

The orbital coordinate system ${\mathrm{O}}_o{\mathrm{X}}_o{\mathrm{Y}}_o{\mathrm{Z}}_o$ of the satellite is shown in Figure 1, and its definition is as follows: the origin O_o is located at the center of mass of the satellite; the O_oZ_o axis points to the center of the earth; O_oY_o points to the negative normal of the orbital plane of the satellite; the O_oX_o, O_oY_o axis and O_oZ_o axis form the right-hand orthogonal coordinate system.

2.1.2. Satellite Body Coordinate System

The body coordinate system of the satellite is represented as ${\mathrm{O}}_{\mathrm{b}}{\mathrm{X}}_{\mathrm{b}}{\mathrm{Y}}_{\mathrm{b}}{\mathrm{Z}}_{\mathrm{b}}$, and the coordinate origin ${\mathrm{O}}_{\mathrm{b}}$ is located at the center of mass of the satellite. When the three-axis attitude of the satellite is zero, the body coordinate system of the satellite coincides with the orbital coordinate system of the satellite.

2.1.3. Nominal Coordinate System

The layout scheme of EP with a mechanical arm adopted in this paper is shown in Figure 2. After the mechanical arm is deployed, the mechanical arm and the solar wing are located in the south and north directions of the GEO satellite, respectively. For the convenience of discussion, it is assumed that during the stable operation of the satellite in orbit, the electric thrust points to the ${\mathrm{Y}}_{\mathrm{b}}$ direction of the satellite.

This kind of EP with a mechanical arm can point to the positive or negative normal direction of the orbital coordinate system, which can achieve efficient NSSK [6]. The thrust vector is slightly deflected to generate tangential thrust during NSSK, and the mean longitude keeping can be achieved at the same time [7]. The thrust vector is slightly deflected to generate radial thrust during NSSK, and the eccentricity keeping can be achieved at the same time [8]. On the basis of the SK function, the position of the end of the mechanical arm is properly deflected, and a wide range of angular momentum unloading can be realized. The motion range of the mechanical arm is small, and the electric thrust only needs to be turned on once per orbit.

If there is only propulsion to the south or north, the incremental speed used for mean longitude keeping will cause the eccentricity vector to diverge in one direction [6,7]. Therefore, the direction of propulsion needs to be alternating north–south during NSSK to prevent the eccentricity vector from diverging in one direction. Definition ${\mathrm{S}}_{\mathsf{\psi }}$ represents the flight state of the satellite. When the satellite is flying forward (the yaw attitude of the satellite $\psi =0°$), ${\mathrm{S}}_{\mathsf{\psi }}>0$, the EP is located on the south side of the satellite, producing a speed increment from south to north, providing the positive normal speed increment of the orbital plane. When the satellite is flying backward (satellite yaw attitude $\psi =180°$), ${\mathrm{S}}_{\mathsf{\psi }}<0$, the electric push is located on the north side of the satellite, producing a speed increment from north to south, providing a negative normal speed increment of the orbital plane.

Due to the limited range of activity of the mechanical arm, it is necessary to change the direction of the mechanical arm with the electric thruster to the south or north by switching the forward flight and backward flight state of the satellite, so as to change the direction of thrust. Therefore, the concept of “nominal coordinate system” is introduced in order to describe the rotation process of the manipulator and the control process of attitude disrotation.

When the satellite is flying forward, the nominal coordinate system coincides with the orbital coordinate system of the satellite. When a satellite is flying backward, the nominal coordinate system is obtained by rotating the orbital coordinate system 180° around its yaw axis. Thus, the quaternion from the J2000 coordinate system to the nominal coordinate system is:

$${\mathit{q}}_{\mathrm{is}}={\mathit{q}}_{\mathrm{io}}\otimes {\mathit{q}}_{\mathrm{os}}$$

(1)

where ${\mathit{q}}_{\mathrm{os}}$ is the attitude quaternion from orbital coordinate system to nominal coordinate system, and its expression is as follows:

$${\mathit{q}}_{\mathrm{os}}=\{\begin{array}{cc}\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]& {\mathrm{S}}_{\mathsf{\psi }}>0\\ \left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right]& {\mathrm{S}}_{\mathsf{\psi }}<0\end{array}$$

(2)

In Formula (1), ${\mathit{q}}_{\mathrm{io}}$ is the attitude quaternion from the J2000 coordinate system to the orbital coordinate system, and its expression is as follows:

$${\mathit{q}}_{\mathrm{io}}=\left[\begin{array}{c}\cos \frac{\Omega +u}{2}\cos \frac{i}{2}\\ \cos \frac{\Omega -u}{2}\sin \frac{i}{2}\\ -\sin \frac{\Omega -u}{2}\sin \frac{i}{2}\\ \sin \frac{\Omega +u}{2}\cos \frac{i}{2}\end{array}\right]\otimes \left[\begin{array}{c}0.5\\ -0.5\\ -0.5\\ 0.5\end{array}\right]\approx \left[\begin{array}{c}\cos \frac{\lambda }{2}\\ 0\\ 0\\ \sin \frac{\lambda }{2}\end{array}\right]\otimes \left[\begin{array}{c}0.5\\ -0.5\\ -0.5\\ 0.5\end{array}\right]$$

(3)

where Ω is the right ascension of ascending node satellite, u is the argument of latitude of the satellite, i is the orbital inclination of the satellite, λ is the right ascension of the satellite.

The function of the nominal coordinate system is to increase the attitude offset of the satellite so that an electric thruster can generate positive or negative normal propulsion on the orbital plane for SK as required.

2.1.4. Mechanical Arm Coordinate System

When the angular momentum of the satellite is kept within a reasonable range in any coordinate system, the angular momentum of the satellite can be unloaded. However, due to the uncertainty of propulsion direction used for position keeping each time, it is complicated to describe the direction of EP during angular momentum unloading in fixed coordinates such as the J2000 coordinate system, radial coordinate system, nominal coordinate system, or satellite body coordinate system. In order to facilitate the analysis of the EP direction, the concept of a “mechanical arm coordinate system ${\mathrm{O}}_{\mathrm{b}}{\mathrm{X}}_{\mathrm{m}}{\mathrm{Y}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$” is introduced. It is defined as follows: the origin of the mechanical arm coordinate system is located at the center of mass $O_{\mathrm{b}}$ of the satellite, the position of propulsion action (center of rotation of the two wrist joints of the mechanical arm) is defined as $O_{\mathrm{m}}$, and the ${\mathrm{Y}}_{\mathrm{m}}$ axis of the mechanical arm coordinate system points to the direction of thrust for SK, as shown in Figure 3.

The mechanical arm needs to complete the tasks of SK and angular momentum unloading while eliminating the influence of the three-axis attitude offset of the satellite. The coordinate system of the mechanical arm is a dynamic coordinate system, which is not only affected by the SK strategy, but also related to the offset attitude of the satellite. The attitude quaternion ${\mathit{q}}_{\mathrm{im}}$ from the J2000 coordinate system to the mechanical arm coordinate system is defined as:

$${\mathit{q}}_{\mathrm{im}}={\mathit{q}}_{\mathrm{is}}\otimes {\mathit{q}}_{\mathrm{sm}}$$

(4)

where ${\mathit{q}}_{\mathrm{sm}}$ is the attitude quaternion from the nominal coordinate system to the mechanical arm coordinate system:

$${\mathit{q}}_{\mathrm{sm}}=\left[\begin{array}{c}\cos \left({\theta }_{\mathrm{t}}/2\right)\\ 0\\ 0\\ \sin \left({\theta }_{\mathrm{t}}/2\right)\end{array}\right]\otimes \left[\begin{array}{c}\cos \left({\theta }_{\mathrm{r}}/2\right)\\ \sin \left({\theta }_{\mathrm{r}}/2\right)\\ 0\\ 0\end{array}\right]$$

(5)

where ${\theta }_t$ is the deflection angle of the EP direction along the ${\mathrm{Z}}_{\mathrm{s}}$ axis of the nominal coordinate system:

$${\theta }_{\mathrm{t}}=\arcsin \frac{\Delta v_{\mathrm{t}}}{\Delta v_{\mathrm{n}}}$$

(6)

When the satellite is flying forward, $\Delta v_{\mathrm{t}}<0$, $\Delta v_{\mathrm{n}}>0$, so ${\theta }_t<0$. This shows that when the satellite is flying forward, the axis ${\mathrm{Y}}_{\mathrm{m}}$ of the mechanical arm coordinate system is obtained by the negative rotation of the ${\mathrm{Y}}_{\mathrm{s}}$ axis of the nominal coordinate system around the ${\mathrm{Z}}_{\mathrm{s}}$ axis, and the mechanical arm is close to the ${\mathrm{X}}_{\mathrm{s}}$ side of the nominal coordinate system. On the contrary, when the satellite is flying backward, $\Delta v_{\mathrm{t}}<0$, $\Delta v_{\mathrm{n}}<0$, so ${\theta }_t<0$. This shows that when the satellite is flying forward, the axis ${\mathrm{Y}}_{\mathrm{m}}$ of the mechanical arm coordinate system is obtained by the positive rotation of the ${\mathrm{Y}}_{\mathrm{s}}$ axis of the nominal coordinate system around the ${\mathrm{Z}}_{\mathrm{s}}$ axis, and the mechanical arm is close to the ${-\mathrm{X}}_{\mathrm{s}}$ side of the nominal coordinate system. In Equation (6), ${\theta }_{\mathit{r}}$ is the deflection angle of the direction of EP along the ${\mathrm{X}}_{\mathrm{s}}$ axis of the nominal coordinate system:

$${\theta }_{\mathrm{r}}=\arcsin \frac{\Delta v_{\mathrm{r}}}{{\mathrm{S}}_{\mathsf{\psi }}\Delta v_{\mathrm{n}}}=\arcsin \frac{\Delta v_{\mathrm{r}}}{\left|\Delta v_{\mathrm{n}}\right|}$$

(7)

The deflection angle of the EP is small, so the three axes of the mechanical arm coordinate system are roughly the same as the three axes of the nominal coordinate system. In Equation (5), both ${\theta }_{\mathrm{r}}$ and ${\theta }_{\mathrm{t}}$ are small angles, and the effect of rotation can be ignored in engineering, so Equation (5) can also be written as:

$${\mathit{q}}_{\mathrm{sm}}=\left[\begin{array}{c}\cos \left({\theta }_{\mathrm{r}}/2\right)\\ \sin \left({\theta }_{\mathrm{r}}/2\right)\\ 0\\ 0\end{array}\right]\otimes \left[\begin{array}{c}\cos \left({\theta }_{\mathrm{t}}/2\right)\\ 0\\ 0\\ \sin \left({\theta }_{\mathrm{t}}/2\right)\end{array}\right]$$

(8)

2.2. Strategy of Unloading Angular Momentum

For GEO satellites with solar panels tracking the sun, the panels are arranged on the north and south sides of the satellite. If the satellite adopts a single-wing sail configuration, under the influence of solar pressure, the angular momentum accumulated in the ${\mathrm{X}}_{\mathrm{b}}{\mathrm{O}}_{\mathrm{b}}{\mathrm{Z}}_{\mathrm{b}}$ plane of the satellite can reach about 40 Nms per day, but the angular momentum accumulated in the satellite ${\pm \mathrm{Y}}_{\mathrm{b}}$ direction is generally less than 1 Nms per day. Therefore, for asymmetric GEO satellites with all-EP, angular momentum unloading is required every day. This means that SK control needs to be performed for a long time each day to achieve the large daily angular momentum offloading requirements.

The semi-major axis of the GEO satellite is about 42,000 km, while the diameter of the Earth is less than 13,000 km. Therefore, the amount of rolling pitch attitude offset required for the satellite to point towards the target on the earth’s surface does not exceed 8.7°. The yaw axis attitude offset required for EWSK does not exceed 5°, and the combined influence of the two on the propulsion direction does not exceed 10°.

In this paper, the angular momentum of the satellite is unloaded by using the propulsion of SK. The unloading of angular momentum is mainly divided into two parts: the first is the large angular momentum unloading in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane of the mechanical arm coordinate system; The second is the small angular momentum unloading in the ${\pm \mathrm{Y}}_{\mathrm{m}}$ direction of the satellite. The method of satellite angular momentum unloading by using propulsion of SK is shown in Figure 4.

Figure 4a is the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane of the mechanical arm coordinate system. This plane will change as the position of the propulsion changes. $O$ is the point of the satellite’s center of mass of projected on the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane (referred to as the origin). ${\mathit{h}}_{xz}^{\mathrm{C}}$ is the angular momentum needed to be unloaded of the satellite in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane of the mechanical arm coordinate system. M is the corresponding propulsion position during angular momentum unloading, and it is also the central position of the EP line movement (referred to as the equivalent position). The direction of OM is the cross of ${\mathit{h}}_{xz}^{\mathrm{C}}$ and ${\mathrm{Y}}_{\mathrm{m}}$. Assuming that the thrust direction is constant along ${\mathrm{Y}}_{\mathrm{m}}$, when the thrust position is M, the unloading amount of angular momentum during the specified SK time is ${\mathit{h}}_{xz}^{\mathrm{C}}$.

As can be seen from Figure 4a, if the propulsion direction is not changed, the angular momentum in the ${\mathrm{Y}}_{\mathrm{m}}$ direction under the mechanical arm coordinate system cannot be unloaded. When ${\mathit{h}}_{xz}^{\mathrm{C}}$ is not affected as much as possible, the angular momentum unloading in the ${\mathrm{Y}}_{\mathrm{m}}$ direction of the mechanical arm coordinate system can be achieved by slightly changing the thrust position and direction at the same time. Specific methods are as follows:

Take any line that goes through M. On this line, two points, A and B, are taken with M as the center, where A is the initial state of the mechanical arm (referred to as initial state A) and B is the final state of the mechanical arm (referred to as final state B). Line MA is the result of line OM rotating ${\beta }_{\mathrm{h}}$ around the ${\mathrm{Y}}_{\mathrm{m}}$ axis. Vectors AC and BD are perpendicular to line AB. Point C is the projection point in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane after the propulsion direction of point A is deflected. Point D is the projection point in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane after the propulsion direction of point B is deflected. Vectors AC and BD are the direction in which the ${\mathrm{Y}}_{\mathrm{m}}$ axis component of propulsion deflected around line AB. Figure 4b shows the three-dimensional space formed before and after the direction change of the thruster. The meaning of the letters in Figure 4b is the same as that in Figure 4a. Vector AC is obtained by ${\mathrm{Y}}_{\mathrm{m}}$ rotating ${\theta }_{\mathrm{m}y}$ around line BA at the starting position A. Vector BD is obtained by ${\mathrm{Y}}_{\mathrm{m}}$ rotating ${\theta }_{\mathrm{m}y}$ around line AB at the final state B.

Assuming that the magnitude F of the propulsion generated by the EP is constant, the length of time that the end of the mechanical arm is in the initial state A and the final state B is the same. In this way, when unloading angular momentum along the ${\mathrm{Y}}_{\mathrm{m}}$ direction, there is almost no disturbed angular momentum in the ${\mathrm{X}}_{\mathrm{m}}$ and ${\mathrm{Z}}_{\mathrm{m}}$ axis. Thus, in order to unload the angular momentum in the Y direction, the position and direction of the end of the mechanical arm are switched between states A and B during the operation of the EP (the position of the EP is switched between A and B, and the direction of the EP is switched between vectors AC and BD).

In the case of the same unloading torque, it is expected that the movement range of the end of the mechanical arm is as small as possible, that is, the minimum value of $\max \left(\left|OA\right|,\left|OB\right|\right)$ is required. $\left|OA\right|$ and $\left|OB\right|$ are, respectively:

$$\begin{array}{l}\left|OA\right|=\sqrt{{\left|AM\right|}^2+{\left|OM\right|}^2+2\left|AM\right|\left|OM\right|\cos {\beta }_{\mathrm{h}}}\\ \left|OB\right|=\sqrt{{\left|BM\right|}^2+{\left|OM\right|}^2-2\left|BM\right|\left|OM\right|\cos {\beta }_{\mathrm{h}}}\end{array}$$

(9)

Then

$$\max \left(\left|OA\right|,\left|OB\right|\right)=\{\begin{array}{cc}\sqrt{{\left|AM\right|}^2+{\left|OM\right|}^2+2\left|AM\right|\left|OM\right|\cos {\beta }_{\mathrm{h}}}& -\pi /2\le {\beta }_{\mathrm{h}}\le \pi /2\\ {\left|AM\right|}^2+{\left|OM\right|}^2-2\left|AM\right|\left|OM\right|\cos {\beta }_{\mathrm{h}}& {\beta }_{\mathrm{h}}>\pi /2,{\beta }_{\mathrm{h}}<-\pi /2\end{array}$$

(10)

According to Equation (11), it can be seen that when ${\beta }_{\mathrm{h}}=\pm \pi /2$, $\max \left(\left|OA\right|,\left|OB\right|\right)$ is the smallest, and the movement range of the end of the mechanical arm is the smallest. Therefore, in subsequent calculations, take ${\beta }_{\mathrm{h}}=\pm \pi /2$.

$R_{OM}$ is the maximum value of vector OM. $R_{OA}$ is the maximum value of vector OA when ${\beta }_{\mathrm{h}}=\pm \pi /2$. It is also the shortest distance between the end of the mechanical arm and the ${\mathrm{Y}}_{\mathrm{m}}$ axis of the mechanical arm coordinate system when the angular momentum is unloaded. Therefore, the maximum arm of force $R_{AM}$ for angular momentum unloading in the ${\pm \mathrm{Y}}_{\mathrm{m}}$ direction is:

$$R_{AM}=\sqrt{{R_{OA}}^2-{R_{OM}}^2}$$

(11)

It can be seen from Equation (11) that when $R_{OA}$ is constant, the longer $R_{OM}$ is, the shorter $R_{AM}$ is. In engineering, appropriate control parameters of angular momentum unloading are selected according to the law of angular momentum variation of satellite.

Let ${\theta }_{\mathrm{m}y\mathrm{MAX}}$ be the maximum deflection angle of EP during angular momentum unloading in the ${\mathrm{Y}}_{\mathrm{m}}$ direction. Then, the ratio of the component $d_{\mathrm{m}y}$ of the arm of force in the ${\mathrm{Y}}_{\mathrm{m}}$ axis generated by angular momentum unloading to the absolute value $\left|{\theta }_{\mathrm{m}y}\right|$ of the deflection angle of the EP is satisfied by Equation (12). It can be seen that the component $d_{\mathrm{m}y}$ of the arm of force in the ${\mathrm{Y}}_{\mathrm{m}}$ axis generated by angular momentum unloading increases linearly as the absolute value $\left|{\theta }_{\mathrm{m}y}\right|$ of the deflection angle of the EP increases.

$$k_{\mathrm{H}y}\triangleq \frac{d_{\mathrm{m}y}}{\left|{\theta }_{\mathrm{m}y}\right|}=\frac{R_{AM}}{{\theta }_{\mathrm{m}y\mathrm{MAX}}}$$

(12)

However, when unloading angular momentum in the ${\mathrm{Y}}_{\mathrm{m}}$ direction, the change in propulsion direction will cause the perturbed moment of force in the AB direction, as shown in Figure 5.

The definition of AB in Figure 5 is consistent with that in Figure 4. In Figure 5a, the thrust vector of the EP is at the initial state A. The effective arm of force corresponding to the unloading moment of force ${\mathit{T}}_{1A}$ generated by force F in the ${\mathrm{Y}}_{\mathrm{m}}$ direction is ${\mathit{r}}_{1A}$. The perturbed arm of force corresponding to the unloading moment of force ${\mathit{T}}_{2A}$ generated by force F in the AB direction is ${\mathit{r}}_{2A}$. In Figure 5b, the thrust vector of the EP is at the terminal state B. The effective arm of force corresponding to the unloading moment of force ${\mathit{T}}_{1\mathrm{B}}$ generated by force F in the ${\mathrm{Y}}_{\mathrm{m}}$ direction is ${\mathit{r}}_{1\mathrm{B}}$. The perturbed arm of force corresponding to the unloading moment of force ${\mathit{T}}_{2B}$ generated by force F in the AB direction is ${\mathit{r}}_{2B}$.

In the two states A and B, the unloading moment of force ${\mathit{T}}_{1A}$ and ${\mathit{T}}_{1\mathrm{B}}$ at A and B are in the same direction and superimposed on each other, respectively. ${\mathit{T}}_{2A}$ and ${\mathit{T}}_{2B}$ are equal in size and opposite in direction, and the resulting perturbed moment of force can cancel each other out in a complete unloading process. Therefore, the unloading moment of force only exists in the ${\mathrm{Y}}_{\mathrm{m}}$ direction. It is worth noting that for the condition ${\mathit{T}}_{2A}$ = ${\mathit{T}}_{2B}$ to hold, the component of the distance in the ${\mathrm{Y}}_{\mathrm{m}}$ direction between the EP and the satellite center of mass must remain constant during the motion of the mechanical arm in the direction AB.

Since the perturbed moment of force ${\mathit{T}}_{2A}$ and ${\mathit{T}}_{2B}$ are larger than the angular momentum unloading moment of force of the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane, the angular momentum of the satellite in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane may be out of tolerance during the unloading process. The solution to this problem is to reduce the single accumulation time of the perturbed moment of force of axis A and B by increasing the switching times $n_{\mathrm{H}y}$ ($n_{\mathrm{H}y}\ge 1$) between the A, B states of the mechanical arm, and then reduce the contribution of the perturbed moment of force A and B to the angular momentum of the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane. The greater the unloading angular momentum required in the ${\mathrm{Y}}_{\mathrm{m}}$ direction, the greater the corresponding $n_{\mathrm{H}y}$. However, the perturbed moment of force corresponding to states A and B are opposite in direction but unequal in magnitude, which will affect the angular momentum unloading effect of the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane and need to be eliminated by the next angular momentum unloading.

2.3. Position and Attitude Calculation of EP with a Mechanical Arm

2.3.1. Connection between EP and a Mechanical Arm

The thrust generated by the EP can be equivalent to the direction of the thrust. The thrust direction passes through the center of the EP, and it can be considered that the equivalent installation position of the EP is on the line of the thrust direction, as shown in Figure 6.

The installation position of the EP is the intersection of the following two lines: ① Thrust vector generated by EP; ② The vertical line between thrust and the rotation axis of the end joint of the robotic arm. If the direction of the EP intersects the rotation axis of the end joint of the robotic arm, the intersection point is the equivalent installation position of the EP. Therefore, the position and attitude control of the end of the robotic arm are decoupled. In addition, the control of the robotic arm has nothing to do with the innovation of this paper, so the relevant content is not explained.

2.3.2. Position and Attitude Calculation of EP in the Mechanical Arm Coordinate System

In Figure 4, the coordinates $\left(x_{\mathrm{m}}^{\mathrm{M}},z_{\mathrm{m}}^{\mathrm{M}}\right)$ of point M in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane are

$$\begin{array}{l}{x}_{\mathrm{m}}^{\mathrm{M}}=-\frac{{\mathit{h}}_{\mathrm{m}z}^{\mathrm{C}}}{F{t}_{\mathrm{NS}}}\\ z_{\mathrm{m}}^{\mathrm{M}}=\frac{{\mathit{h}}_{\mathrm{m}x}^{\mathrm{C}}}{F{t}_{\mathrm{NS}}}\end{array}$$

(13)

where $h_{\mathrm{m}x}^{\mathrm{C}}$ and $h_{\mathrm{m}z}^{\mathrm{C}}$ are the components of ${\mathit{h}}_{xz}^{\mathrm{C}}$ on the ${\mathrm{X}}_{\mathrm{m}}$ axis and ${\mathrm{Z}}_{\mathrm{m}}$ axis of the mechanical arm coordinate system, respectively; $t_{\mathrm{NS}}$ is the duration of orbital inclination keeping using EP. The formula for ${\mathit{h}}_{xz}^{\mathrm{C}}$ is

$${\mathit{h}}_{xz}^{\mathrm{C}}\left[\begin{array}{c}{h}_{\mathrm{m}x}^{\mathrm{C}}\\ h_{\mathrm{m}y}^{\mathrm{C}}\\ h_{\mathrm{m}z}^{\mathrm{C}}\end{array}\right]={\mathit{A}}_{\mathrm{m}\mathrm{i}}\left[\begin{array}{c}{h}_{\mathrm{i}x}^{\mathrm{C}}\\ h_{\mathrm{i}y}^{\mathrm{C}}\\ h_{\mathrm{i}z}^{\mathrm{C}}\end{array}\right]$$

(14)

where ${\mathit{A}}_{\mathrm{m}\mathrm{i}}$ is the attitude transformation matrix from the J2000 coordinate system to the mechanical arm coordinate system. ${\left[\begin{array}{ccc}{h}_{\mathrm{i}x}^{\mathrm{C}}& h_{\mathrm{i}y}^{\mathrm{C}}& h_{\mathrm{i}z}^{\mathrm{C}}\end{array}\right]}^{\mathrm{T}}$ is the angular momentum that needs to be unloaded of the satellite in the J2000 coordinate system. Its calculation formula is

$$\left[\begin{array}{c}{h}_{\mathrm{i}x}^{\mathrm{C}}\\ h_{\mathrm{i}y}^{\mathrm{C}}\\ h_{\mathrm{i}z}^{\mathrm{C}}\end{array}\right]=\left[\begin{array}{c}-0.5T{\dot{h}}_{\mathrm{i}x}-H_{\mathrm{i}x}\\ -0.5T{\dot{h}}_{\mathrm{i}y}-H_{\mathrm{i}y}\\ -0.5T{\dot{h}}_{\mathrm{i}z}-H_{\mathrm{i}z}\end{array}\right]$$

(15)

where ${\left[\begin{array}{ccc}{\dot{h}}_{\mathrm{i}x}& {\dot{h}}_{\mathrm{i}y}& {\dot{h}}_{\mathrm{i}z}\end{array}\right]}^{\mathrm{T}}$ is the perturbed moment of force of the satellite in the J2000 coordinate system.

In order to make the angular momentum before and after unloading evenly distributed in center zero as far as possible, it is necessary to calculate the angular momentum to be unloaded in the way of “negative offset”. As shown in Figure 7, Point A is the angular momentum of the satellite in the J2000 coordinate system when the angular momentum unloading strategy is generated. Point C is the angular momentum of the satellite in the J2000 coordinate system under the action of perturbed moment of force of solar light pressure. O is the desired angular momentum envelope center. Point D is the angular momentum that needs to be unloaded this time. AC and DO are the direction of solar pressure perturbed moment of force. When the angular momentum is unloaded to point D, under the action of perturbed moment of force of the solar pressure, the angular momentum of the satellite will pass through point O, and the trajectory of the angular momentum will be evenly distributed at both ends of point O.

Assuming that the satellite is equipped with a single-wing solar panel, because the center of solar pressure of the satellite does not coincide with the center of mass of the satellite, the solar panel of the satellite continues to point to the sun, which will produce approximately constant perturbed moment of force under the action of solar light pressure. Since the sun moves periodically in the J2000 coordinate system, the perturbed moment of force also changes periodically. Therefore, ${\left[\begin{array}{ccc}{\dot{h}}_{\mathrm{i}x}& {\dot{h}}_{\mathrm{i}y}& {\dot{h}}_{\mathrm{i}z}\end{array}\right]}^{\mathrm{T}}$ in Equation (15) can be regarded as a constant value in a short period of time, and its calculation is as follows:

$$\left[\begin{array}{c}{\dot{h}}_{\mathrm{i}x}\\ {\dot{h}}_{\mathrm{i}y}\\ {\dot{h}}_{\mathrm{i}z}\end{array}\right]=\{\begin{array}{cc}{\mathit{T}}_{\mathrm{P}}& {\mathrm{S}}_{\mathsf{\psi }}>0\\ {{\mathit{T}}_{\mathrm{P}}}^{\prime }& {\mathrm{S}}_{\mathsf{\psi }}<0\end{array}$$

(16)

where ${\mathit{T}}_{\mathrm{P}}$ is the perturbed moment of force of solar light pressure in the J2000 coordinate system when the satellite is flying forward. ${{\mathit{T}}_{\mathrm{P}}}^{\prime }$ is the perturbed moment of force of solar light pressure in the J2000 coordinate system when the satellite is flying backward, and ${{\mathit{T}}_{\mathrm{P}}}^{\prime }=-{\mathit{T}}_{\mathrm{P}}$. In Equation (15), ${\left[\begin{array}{ccc}{H}_{\mathrm{i}x}& H_{\mathrm{i}y}& H_{\mathrm{i}z}\end{array}\right]}^{\mathrm{T}}$ is the angular momentum of the satellite in the J2000 coordinate system. Its calculation formula is as follows:

$$\left[\begin{array}{c}{H}_{\mathrm{i}x}\\ H_{\mathrm{i}y}\\ H_{\mathrm{i}z}\end{array}\right]={\mathit{A}}_{\mathrm{ib}}\left[\begin{array}{c}{H}_{\mathrm{b}x0}\\ H_{\mathrm{b}y0}\\ H_{\mathrm{b}z0}\end{array}\right]+\left[\begin{array}{c}{\dot{h}}_{\mathrm{i}x}\\ {\dot{h}}_{\mathrm{i}y}\\ {\dot{h}}_{\mathrm{i}z}\end{array}\right]\Delta t_{\mathrm{h}}$$

(17)

where ${\mathit{A}}_{\mathrm{ib}}$ is the attitude transformation matrix from the satellite body coordinate system to the J2000 coordinate system. ${\left[\begin{array}{ccc}{H}_{\mathrm{b}x0}& H_{\mathrm{b}y0}& H_{\mathrm{b}z0}\end{array}\right]}^{\mathrm{T}}$ is the angular momentum of the satellite in the satellite body coordinate system when the angular momentum unloading strategy is generated. $\Delta t_{\mathrm{h}}$ is the time difference between the time when the angular momentum unloading strategy is generated and the intermediate time when the next SK pulse is applied. Its calculation formula is

$$\Delta t_{\mathrm{h}}=\frac{\mathrm{mod}\left({\lambda }_{\mathrm{tm}}-\lambda ,2\pi \right)}{n}$$

(18)

where ${\lambda }_{\mathrm{tm}}$ is the orbital right ascension corresponding to the midpoint of EP operation. $\lambda$ is the orbital right ascension at the time when the angular momentum unloading strategy is generated. n is the average orbital angular velocity.

Considering the limited range of motion of the end of the mechanical arm, Equation (13) can be written as:

$$\begin{array}{l}{x}_{\mathrm{m}}^M=\frac{-h_{\mathrm{m}z}^{\mathrm{C}}\min \left(\sqrt{{\left(h_{\mathrm{m}z}^{\mathrm{C}}\right)}^2+{\left(h_{\mathrm{m}x}^{\mathrm{C}}\right)}^2},F{t}_{\mathrm{N}\mathrm{S}}{R}_{OM}\right)}{F{t}_{\mathrm{NS}}\sqrt{{\left(h_{\mathrm{m}z}^{\mathrm{C}}\right)}^2+{\left(h_{\mathrm{m}x}^{\mathrm{C}}\right)}^2}}\\ z_{\mathrm{m}}^M=\frac{h_{\mathrm{m}x}^{\mathrm{C}}\min \left(\sqrt{{\left(h_{\mathrm{m}z}^{\mathrm{C}}\right)}^2+{\left(h_{\mathrm{m}x}^{\mathrm{C}}\right)}^2},F{t}_{\mathrm{N}\mathrm{S}}{R}_{OM}\right)}{F{t}_{\mathrm{NS}}\sqrt{{\left(h_{\mathrm{m}z}^{\mathrm{C}}\right)}^2+{\left(h_{\mathrm{m}x}^{\mathrm{C}}\right)}^2}}\end{array}$$

(19)

The electric thrust deflection angle ${\theta }_{\mathrm{m}y}$ is calculated as follows:

$${\theta }_{\mathrm{m}y}=\mathrm{sgn}\left(h_{\mathrm{m}y}^{\mathrm{C}}\right)\sqrt{\frac{h_{\mathrm{m}y}^{\mathrm{C}}}{k_{\mathrm{H}y}F{t}_{\mathrm{NS}}}}$$

(20)

where $h_{\mathrm{m}y}^{\mathrm{C}}$ is the unloading angular momentum in the ${\mathrm{Y}}_{\mathrm{m}}$ direction.

The value range of ${\theta }_{\mathrm{m}y}$ is ${\theta }_{\mathrm{m}y}\in \left(-90°,90°\right)$. The polarity of ${\theta }_{\mathrm{m}y}$ indicates the direction of EP. For example, when ${\theta }_{\mathrm{m}y}>0$, it means that the EP rotates positively along the MA direction at point A, and the EP rotates positively along the MB direction at point B. In specific applications, ${\theta }_{\mathrm{m}y}$ should be limited to an appropriate range according to actual engineering constraints. Formula (20) can be written as:

$${\theta }_{\mathrm{m}y}=\{\begin{array}{cc}\mathrm{sgn}\left(h_{\mathrm{m}y}^{\mathrm{C}}\right)\sqrt{h_{\mathrm{m}y}^{\mathrm{C}}/\left(k_{\mathrm{H}y}F{t}_{\mathrm{NS}}\right)}& \sqrt{h_{\mathrm{m}y}^{\mathrm{C}}/\left(k_{\mathrm{H}y}F{t}_{\mathrm{NS}}\right)}\le {\theta }_{\mathrm{myMAX}}\\ \mathrm{sgn}\left(h_{\mathrm{m}y}^{\mathrm{C}}\right){\theta }_{\mathrm{m}y\mathrm{MAX}}& \sqrt{h_{\mathrm{m}y}^{\mathrm{C}}/\left(k_{\mathrm{H}y}F{t}_{\mathrm{NS}}\right)}>{\theta }_{\mathrm{myMAX}}\end{array}$$

(21)

${\theta }_{\mathrm{myMAX}}$ is the maximum value of the electric deflection angle when the angular momentum in the direction ${\mathrm{Y}}_{\mathrm{m}}$ is unloaded, which is positive.

In Figure 4, the coordinates $\left(x_{\mathrm{m}}^A,z_{\mathrm{m}}^A\right)$ of point A in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane in the mechanical arm coordinate system are

$$\begin{array}{l}{x}_{\mathrm{m}}^A=x_{\mathrm{m}}^M-k_{\mathrm{H}y}\left|{\theta }_{\mathrm{m}y}\right|\sin {\beta }_{\mathrm{h}}\\ z_{\mathrm{m}}^A=z_{\mathrm{m}}^M-k_{\mathrm{H}y}\left|{\theta }_{\mathrm{m}y}\right|\cos {\beta }_{\mathrm{h}}\end{array}$$

(22)

The coordinates $\left(x_{\mathrm{m}}^B,z_{\mathrm{m}}^B\right)$ of point B in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane in the mechanical arm coordinate system are

$$\begin{array}{l}{x}_{\mathrm{m}}^B=x_{\mathrm{m}}^M+k_{\mathrm{H}y}\left|{\theta }_{\mathrm{m}y}\right|\sin {\beta }_{\mathrm{h}}\\ z_{\mathrm{m}}^B=z_{\mathrm{m}}^M+k_{\mathrm{H}y}\left|{\theta }_{\mathrm{m}y}\right|\cos {\beta }_{\mathrm{h}}\end{array}$$

(23)

The unit vector ${\mathit{V}}_{\mathrm{m}}^A$ of A in the mechanical arm system is

$${\mathit{V}}_{\mathrm{m}}^A=\cos {\theta }_{\mathrm{m}y}{\mathit{V}}_{\mathrm{m}}^0+\left(1-\cos {\theta }_{\mathrm{m}y}\right)\left({\mathit{V}}_{\mathrm{m}}^{MA}\cdot {\mathit{V}}_{\mathrm{m}}^0\right){\mathit{V}}_{\mathrm{m}}^{MA}+\sin {\theta }_{\mathrm{m}\mathrm{y}}\left({\mathit{V}}_{\mathrm{m}}^{MA}\times {\mathit{V}}_{\mathrm{m}}^0\right)$$

(24)

where ${\mathit{V}}_m^0={\left(\begin{array}{ccc}0& 1& 0\end{array}\right)}^{\mathrm{T}}$ is the initial direction of the EP in the mechanical arm coordinate system. The expression for ${\mathit{V}}_{\mathrm{m}}^{MA}$ is:

$${\mathit{V}}_{\mathrm{m}}^{MA}=\left[\begin{array}{c}\sin \left(\frac{\pi }{2}+{\beta }_{\mathrm{h}0}\right)\\ 0\\ \cos \left(\frac{\pi }{2}+{\beta }_{\mathrm{h}0}\right)\end{array}\right]$$

(25)

where ${\beta }_{\mathrm{h}0}$ is the angle between ${\mathrm{OZ}}_{\mathrm{m}}$ and OM in the ${\mathrm{X}}_{\mathrm{m}}{\mathrm{O}}_{\mathrm{m}}{\mathrm{Z}}_{\mathrm{m}}$ plane, expressed as:

$${\beta }_{\mathrm{h}0}=\{\begin{array}{cc}\arctan \left(h_{\mathrm{m}x}^{\mathrm{C}}/h_{\mathrm{m}z}^{\mathrm{C}}\right)& h_{\mathrm{m}z}^{\mathrm{C}}>0\\ \pi +\arctan \left(h_{\mathrm{m}x}^{\mathrm{C}}/h_{\mathrm{m}z}^{\mathrm{C}}\right)& h_{\mathrm{m}z}^{\mathrm{C}}<0\\ 0& h_{\mathrm{m}z}^{\mathrm{C}}=0\end{array}$$

(26)

The unit vector ${\mathit{V}}_{\mathrm{m}}^B$ of B in the mechanical arm system is

$${\mathit{V}}_{\mathrm{m}}^B=\cos {\theta }_{\mathrm{m}y}{\mathit{V}}_{\mathrm{m}}^0+\left(1-\cos {\theta }_{\mathrm{m}y}\right)\left({\mathit{V}}_{\mathrm{m}}^{MB}\cdot {\mathit{V}}_{\mathrm{m}}^0\right){\mathit{V}}_{\mathrm{m}}^{MB}+\sin {\theta }_{\mathrm{m}y}\left({\mathit{V}}_{\mathrm{m}}^{MB}\times {\mathit{V}}_{\mathrm{m}}^0\right)$$

(27)

where

$${\mathit{V}}_{\mathrm{m}}^{MB}=-{\mathit{V}}_{\mathrm{m}}^{MA}$$

(28)

2.3.3. Position and Attitude Calculation of EP in the Satellite Body Coordinate System

In order to facilitate the simulation and evaluation of the motion range of the mechanical arm in the satellite body coordinate system, the position and attitude of the EP in the satellite body coordinate system should also be calculated. The coordinates $\left(x_b^A,z_b^A\right)$ of point A in the ${\mathrm{X}}_{\mathrm{b}}{\mathrm{O}}_{\mathrm{b}}{\mathrm{Z}}_{\mathrm{b}}$ plane in the satellite body coordinate system are

$$\left(\begin{array}{c}{x}_{\mathrm{b}}^A\\ y_{\mathrm{b}}^A\\ z_{\mathrm{b}}^A\end{array}\right)={\left({\mathit{A}}_{\mathrm{mb}}\right)}^{-1}\left(\begin{array}{c}{x}_{\mathrm{m}}^A\\ 0\\ z_{\mathrm{m}}^A\end{array}\right)$$

(29)

The coordinates $\left(x_b^B,z_b^B\right)$ of point B in the ${\mathrm{X}}_{\mathrm{b}}{\mathrm{O}}_{\mathrm{b}}{\mathrm{Z}}_{\mathrm{b}}$ plane in the satellite body coordinate system are

$$\left(\begin{array}{c}{x}_{\mathrm{b}}^B\\ y_{\mathrm{b}}^B\\ z_{\mathrm{b}}^B\end{array}\right)={\left({\mathit{A}}_{\mathrm{mb}}\right)}^{-1}\left(\begin{array}{c}{x}_{\mathrm{m}}^B\\ 0\\ z_{\mathrm{m}}^B\end{array}\right)$$

(30)

The unit vector ${\mathit{V}}_{\mathrm{b}}^A$ of A in the satellite body coordinate system is:

$${\mathit{V}}_{\mathrm{b}}^A={\left({\mathit{A}}_{\mathrm{mb}}\right)}^{-1}{\mathit{V}}_{\mathrm{m}}^A$$

(31)

The unit vector ${\mathit{V}}_{\mathrm{b}}^B$ of B in the satellite body coordinate system is:

$${\mathit{V}}_{\mathrm{b}}^B={\left({\mathit{A}}_{\mathrm{mb}}\right)}^{-1}{\mathit{V}}_{\mathrm{m}}^B$$

(32)

In Equations (29)–(32), ${\mathit{A}}_{\mathrm{mb}}$ is the attitude transfer matrix from the mechanical arm coordinate system to the satellite body coordinate system.

2.4. Angular Momentum Unloading in the Case of Attitude Offset

In general, when the satellite has a three-axis attitude offset (attitude maneuver can also be regarded as attitude offset), the direction of thrust used for SK will also be offset. In this paper, a mechanical arm carrying electric thrust is used to control the thrust direction as expected during the three-axis attitude offset.

Take the GEO satellite’s Z axis pointing to the center of the earth as an example. The diagrams of a satellite pointing at a target on the earth’s surface as it is flying forward and backward are shown in Figure 8 and Figure 9, respectively.

Assume that the earth is a standard spherical sphere (without considering the oblate of the earth), O is the center of the earth, B is the object to be observed on the earth’s surface, $O_{\mathrm{b}}$ is the center of mass of the satellite, ${\mathrm{X}}_{\mathrm{s}}{,\mathrm{Y}}_{\mathrm{s}}{,\mathrm{Z}}_{\mathrm{s}}$ are, respectively, the three-axis direction of the nominal coordinate system, $\varphi$ is the latitude of the object to be observed on the earth’s surface, and $\Delta \lambda$ is the relative longitude between the object to be observed on the earth’s surface and the satellite’s subsatellite point.

When the ${\mathrm{Z}}_{\mathrm{b}}$ axis of the satellite body coordinate system points to the ground target B, its three-axis coordinate system can be regarded as the satellite rotating $\gamma$ from the nominal coordinate system around the space line $O_{\mathrm{b}}A$. The direction of line $O_{\mathrm{b}}A$ is:

$$O_{\mathrm{b}}A=O_{\mathrm{b}}O\times OB$$

(33)

In the satellite body coordinate system, the expression of $O_{\mathrm{b}}O$ is:

$$O_{\mathrm{b}}O=\left[\begin{array}{c}0\\ 0\\ a\end{array}\right]$$

(34)

In the satellite body coordinate system, the expression of $OB$ is:

$$OB=\left[\begin{array}{c}{r}_{\mathrm{e}}\cos \varphi \sin \Delta \lambda \\ -r_{\mathrm{e}}\sin \varphi \\ -r_{\mathrm{e}}\cos \varphi \cos \Delta \lambda \end{array}\right]$$

(35)

Substitute $O_{\mathrm{b}}O$ and $OB$ into Formula (9) to obtain

$$O_{\mathrm{b}}A=\left[\begin{array}{c}a{r}_{\mathrm{e}}\sin \varphi \\ a{r}_{\mathrm{e}}\cos \varphi \sin \Delta \lambda \\ 0\end{array}\right]$$

(36)

In the satellite body coordinate system, the expression of $O_{\mathrm{b}}B$ is:

$$O_{\mathrm{b}}B=\left[\begin{array}{c}{r}_{\mathrm{e}}\cos \varphi \sin \Delta \lambda \\ -r_{\mathrm{e}}\sin \varphi \\ a-r_{\mathrm{e}}\cos \varphi \cos \Delta \lambda \end{array}\right]$$

(37)

According to the law of cosine, $\gamma$ can be calculated as:

$$\gamma =\arccos \frac{{\left|O_{\mathrm{b}}O\right|}^2+{\left|O_{\mathrm{b}}B\right|}^2-{\left|OB\right|}^2}{2\left|O_{\mathrm{b}}O\right|\left|O_{\mathrm{b}}B\right|}$$

(38)

The attitude quaternion ${\mathit{q}}_{\mathrm{sb}}$ from the nominal coordinate system to the satellite body coordinate system is

$${\mathit{q}}_{\mathrm{sb}}=\left[\begin{array}{c}\cos \frac{\gamma }{2}\\ \frac{\sin \varphi }{\sqrt{{\left(\sin \varphi \right)}^2+{\left(\cos \varphi \sin \Delta \lambda \right)}^2}}\sin \frac{\gamma }{2}\\ \frac{\cos \varphi \sin \Delta \lambda }{\sqrt{{\left(\sin \varphi \right)}^2+{\left(\cos \varphi \sin \Delta \lambda \right)}^2}}\sin \frac{\gamma }{2}\\ 0\end{array}\right]$$

(39)

By combining Equations (1) and (39), the attitude quaternion ${\mathit{q}}_{\mathrm{ob}}$ from the satellite orbital coordinate system to the body coordinate system can be obtained as follows:

$${\mathit{q}}_{\mathrm{ob}}=\{\begin{array}{cc}\left[\begin{array}{c}\cos \frac{\gamma }{2}\\ \frac{\sin \varphi }{\sqrt{{\left(\sin \varphi \right)}^2+{\left(\cos \varphi \sin \Delta \lambda \right)}^2}}\sin \frac{\gamma }{2}\\ \frac{\cos \varphi \sin \Delta \lambda }{\sqrt{{\left(\sin \varphi \right)}^2+{\left(\cos \varphi \sin \Delta \lambda \right)}^2}}\sin \frac{\gamma }{2}\\ 0\end{array}\right]& {\mathrm{S}}_{\mathsf{\psi }}>0\\ \left[\begin{array}{c}0\\ \frac{\cos \varphi \sin \Delta \lambda }{\sqrt{{\left(\sin \varphi \right)}^2+{\left(\cos \varphi \sin \Delta \lambda \right)}^2}}\sin \frac{\gamma }{2}\\ -\frac{\sin \varphi }{\sqrt{{\left(\sin \varphi \right)}^2+{\left(\cos \varphi \sin \Delta \lambda \right)}^2}}\sin \frac{\gamma }{2}\\ \cos \frac{\gamma }{2}\end{array}\right]& {\mathrm{S}}_{\mathsf{\psi }}<0\end{array}$$

(40)

In the nominal coordinate system, the three-axis attitude of the satellite is the corresponding nominal forward or backward flight attitude. Therefore, angular momentum unloading can be accomplished by executing the SK strategy [16,17,18] in the nominal coordinate system. The conversion matrix from the nominal coordinate system to the satellite body coordinate system is the conversion matrix of the satellite attitude offset. Therefore, the attitude quaternion ${\mathit{q}}_{\mathrm{bm}}$ from the satellite body coordinate system to the mechanical arm coordinate system is:

$${\mathit{q}}_{\mathrm{bm}}={\mathit{q}}_{\mathrm{bo}}\otimes {\mathit{q}}_{\mathrm{om}}={{\mathit{q}}_{\mathrm{ob}}}^{-1}\otimes {\mathit{q}}_{\mathrm{om}}$$

(41)

where ${\mathit{q}}_{\mathrm{om}}$ is the attitude quaternion from the orbital coordinate system to the mechanical arm coordinate system:

$${\mathit{q}}_{\mathrm{om}}={{\mathit{q}}_{\mathrm{io}}}^{-1}\otimes {\mathit{q}}_{\mathrm{im}}$$

(42)

In addition, ${\mathit{q}}_{\mathrm{bm}}$ can also be used to convert the angular momentum of the satellite body coordinate system to the mechanical arm coordinate system and then unload the satellite angular momentum under the mechanical arm coordinate system.

The above method can be applied to any satellite attitude within the reach range of the mechanical arm. The attitude maneuver can be regarded as a special attitude offset (that is, the attitude angle of the satellite is continuously offset), so the method in this paper is also applicable to the angular momentum unloading during the attitude maneuver.

2.5. Fuel Consumption Analysis for Angular Momentum Unloading

The angular momentum unloading strategy adopted in this paper is as follows: in an SK arc, the position and attitude of the EP are changed by the mechanical arm, and the unloading torque in the ${\mathrm{Y}}_{\mathrm{m}}$ direction is generated by the SK, so as to realize the angular momentum unloading of the satellite. Therefore, the deflection angle ${\theta }_{my}$ of the EP determines the efficiency of SK. The smaller the absolute value of the deflection angle, the smaller the loss of thrust for SK.

Take $\left|{\theta }_{my}\right|=5^{\circ }$ as an example, then the efficiency of thrust is $\cos {\theta }_{my}=0.9962$. As a result, the loss of SK efficiency due to angular momentum unloading is at most 0.4%. According to the literature [16,17,18], the annual speed increment required for SK does not exceed 50 m/s. As a result, the required velocity increment for angular momentum unloading is less than 0.2 m/s per year.

3. Simulation Result

According to the analysis in Section 4, in order to completely cancel the perturbed moment of force generated in states A and B in Figure 5, the perturbed moment of force at A and B are required to be equal in magnitude. This requires that the distance between the mechanical arm at point A $y_{\mathrm{b}}^A$ and the center of mass of the satellite and the distance between the mechanical arm at point B $y_{\mathrm{b}}^B$ and the center of mass of the satellite are equal during the process of moving the end of the mechanical arm. The larger $y_{\mathrm{b}}^A$ and $y_{\mathrm{b}}^B$ are, the larger the perturbed moment of force generated during angular momentum unloading. If $y_{\mathrm{b}}^A$ and $y_{\mathrm{b}}^B$ are too large, the angular momentum of the satellite is out of tolerance. Therefore, in the simulation, it can be assumed that both $y_{\mathrm{b}}^A$ and $y_{\mathrm{b}}^B$ are 4 m according to satellite configuration constraints. Suppose that the accumulated angular momentum in the ${\mathrm{Y}}_{\mathrm{b}}$ direction of the satellite body coordinate system is 2 Nms a day. Then, the angular momentum accumulated in the ${\mathrm{X}}_{\mathrm{b}}{\mathrm{O}}_{\mathrm{b}}{\mathrm{Z}}_{\mathrm{b}}$ plane of the satellite body coordinate system is about 46 Nms per day. Let ${\theta }_{\mathrm{myMAX}}=5^{\circ }$, $n_{\mathrm{H}y}=3$, and other initial conditions for simulation are shown in

Table 1. Orbital parameters used in simulation.

Name	Parameter
Coordinate System	J2000 Coordinate System
Orbital Epoch $T_0$ (UTC)	1 August 2025 12:00:00
Semi-major Axis $a$ (km)	42,166.3
Eccentricity $e$	0.0001
Inclination $i$ (°)	0.08
Argument of Perigee $\omega$ (°)	0
RAAN $\Omega$ (°)	359.989
Mean Anomaly $M$ (°)	251.361
Mean Longitude $l$ (°)	121.019
Mass of Satellite $m$ (kg)	3000
$\mathrm{Solar} \mathrm{Luminous} \mathrm{Pressure} \mathrm{Coefficient} C_R$	1.3

.

3.1. Simulation of Angular Momentum Unloading without Attitude Offset

The simulation results of angular momentum unloading without attitude offset are shown in Figure 10.

As can be seen from Figure 10a,b, the position of EP in the ${\mathrm{X}}_{\mathrm{b}}$ direction is ±0.60 m, and the position in the ${\mathrm{Z}}_{\mathrm{b}}$ direction is ±0.20 m. Assuming that the distance between the EP and the satellite’s center of mass is about 4 m, the EP is active in a square cone with the ${\mathrm{Y}}_{\mathrm{b}}$ axis as the central axis and the east–west semi-cone angle is about 8.5°, and the north–south semi-cone angle is about 2.8°. It can be seen from Figure 10c,d that the maximum value of the angular momentum amplitude of the satellite is 40 Nms, which proves that the strategy of angular momentum unloading proposed in this paper is effective.

3.2. Simulation of Angular Momentum Unloading with Attitude Offset

The simulation results of angular momentum unloading without attitude offset are shown in Figure 11.

As can be seen from Figure 11a,b, the position of EP in the ${\mathrm{X}}_{\mathrm{b}}$ direction is ±0.60 m, and the position in the ${\mathrm{Z}}_{\mathrm{b}}$ direction is ±0.75 m. Assuming that the distance between the electric thrust and the satellite’s center of mass is about 4 m, the electric thrust is active in a square cone with the ${\mathrm{Y}}_{\mathrm{b}}$ axis as the central axis and the east–west semi-cone angle is about 8.5°, and the north–south semi-cone angle is about 10.6°. It can be seen from Figure 11c,d that the maximum value of the angular momentum amplitude of the satellite is 40 Nms, which proves that the strategy of angular momentum unloading proposed in this paper is effective.

4. Discussion

By comparing Figure 10a,b and Figure 11a,b, it can be seen that when the attitude is offset, the position and attitude range of the end of the robotic arm in the ${\mathrm{Z}}_{\mathrm{b}}$ direction of the satellite body coordinate system are significantly larger. This is because when the satellite is required to point to objects near the North and South poles, the attitude of the rolling axis of the satellite is offset, and the angular momentum unloading causes the position and attitude range of the EP in the ${\mathrm{Z}}_{\mathrm{b}}$ direction to become larger.

Due to fuel consumption and other reasons, the actual center of mass of the satellite will be different from the theoretical value. The deviation of the center of mass is generally distributed in a spherical region with a radius of about 0.01 m. When the angular momentum unloading is carried out with an electric push of 80 mN, the maximum thrust action time is 7688 s, and the deviation of angular momentum unloading is about 4.6 Nms. Generally, the angular momentum envelope capacity of large GEO satellites is above 50 Nms, so this deviation of angular momentum unloading can be ignored. The deviation of angular momentum unloading caused by the deviation of the center of mass is mainly manifested as the constant deviation of the center value of the angular momentum of the satellite, which reduces the safety margin of the angular momentum envelope of the satellite. Therefore, it is necessary to correct this deviation. There are two main methods of correction: First, the angular momentum unloading strategy is generated using the precise satellite’s center of mass. This method can completely solve the above problem of angular momentum unloading deviation. Generally, three-axis attitude control by chemical propulsion is used to achieve accurate estimation of the center of mass, but the method is not operable. Second, the target value of this angular momentum unloading can be obtained by using the results of the last angular momentum unloading. This method is relatively simple. However, with the change of seasons, the unloading deviation of angular momentum may have a yearly cycle change, and the target value of angular momentum unloading needs to be updated regularly.

5. Conclusions

This paper presents a method of angular momentum unloading using a mechanical arm with EP. In this method, the satellite can be unloaded with large angular momentum by using an electric thruster in one SK arc. The problems of weak angular momentum unloading capacity and high fuel consumption of all-electric GEO satellites are solved. At the same time, the method can realize angular momentum unloading under the attitude offset of the satellite, which is conducive to the continuous and high-quality operation of the load on board. When the demand for angular momentum unloading along the thrust direction increases, the angle from which the thrust direction deviates from the center direction increases.

The angular momentum unloading algorithm proposed in this paper has a certain universality in the range of angular momentum unloading capacity. This will not only reduce the efficiency of SK but also cause the perturbed moment of force perpendicular to the thrust direction to increase, and the number of switching between the initial state and the final state will increase correspondingly. Therefore, the smaller the angular momentum unloading demand along the thrust direction, the better the adaptability of the angular momentum unloading scheme proposed in this paper.