Next Article in Journal
The Need for Multi-Sensor Data Fusion in Structural Health Monitoring of Composite Aircraft Structures
Next Article in Special Issue
Orbit-Injection Strategy and Trajectory-Planning Method of the Launch Vehicle under Power Failure Conditions
Previous Article in Journal
Icon Design for Representing Safety-Critical Aircraft Functions to Support Supervisory Control of Remotely Piloted Aircraft Systems
Previous Article in Special Issue
A Beam Search-Based Channel Allocation Method for Interference Mitigation of NGSO Satellites with Multi-Beam Antennas
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Geostationary Station-Keeping of Electric-Propulsion Satellite Equipped with Robotic Arms

1
School of Astronomy and Space Science, Nanjing University, Nanjing 210093, China
2
School of Aeronautics and Astronautics, Sun Yat-Sen University, Shenzhen 518107, China
3
China Astronaut Research and Training Center, Beijing 100094, China
*
Author to whom correspondence should be addressed.
Aerospace 2022, 9(4), 182; https://doi.org/10.3390/aerospace9040182
Submission received: 21 February 2022 / Revised: 22 March 2022 / Accepted: 23 March 2022 / Published: 28 March 2022
(This article belongs to the Special Issue Recent Advances in Spacecraft Dynamics and Control)

Abstract

:
We propose two approaches based on feedforward control and model-predictive control, respectively, to solve the station-keeping problem of an electric-propulsion geostationary Earth orbit (GEO) satellite, whose thrusters are mounted on two robotic arms on its anti-nadir face. This novel configuration enables a wider range of thrust direction, making it possible to regard the thrust direction as control variables. To solve this control problem, we present the quick feedforward controller (QFFC) and the fuel-optimal model predictive controller (FOMPC). The QFFC is developed based on the analysis of GEO dynamics and the thruster configuration. The FOMPC applies an optimization algorithm to solve the nonlinear model predictive control (NLMPC) problem with the initial value given by the QFFC. Numerical simulations suggest that both controllers could achieve stable station-keeping over multiple objective elements with fewer thrusters and fewer maneuvers. The QFFC has higher control accuracy and lower computational requirements than the FOMPC, whereas the FOMPC could significantly save fuel consumption. The robustness assessment and other discussions of the controllers are also presented.

1. Introduction

Geostationary Earth orbit (GEO) satellites have been playing an important role in applications such as communications, Earth environment, navigation, and aircraft control for decades [1]. To help adjacent satellites avoid colliding with each other, the number of satellites that could be deployed safely around a certain longitude is limited [2]. High-accuracy station-keeping control could alleviate this problem by reducing the required minimum distance or enabling formation flying [3].
Station-keeping by using chemical propulsion is a mature technology and has been applied in plenty of missions [4,5], in which the ground-control center generates the control command by measuring the deviation of the satellite. Due to the nature of the chemical propulsion, chemical-propulsion GEO satellites execute station-keeping maneuvers every two weeks conventionally, resulting in a significant decline in the station-keeping accuracy. In order to solve the problems caused by chemical propulsion, the studies on new types of propulsion, especially electric propulsion, can be traced back to the 1950s [6]. Today, the electric-propulsion satellite has been capable of completing various kinds of space missions [7,8,9], making it possible to form an all-electric satellite platform [10]. The superiority of electric propulsion consists in the high specific impulse, which could significantly reduce the fuel consumption over a lifetime of 12–15 years of GEO satellites [11]. In the past 30 years, multiple electric propulsion systems, such as BPT-4000 [12], SPT-100 [13], and PPS-1350 [14], have been successfully deployed in missions, and the Boeing company also introduced their world’s first all-electric satellite design 702SP in 2012 [10]. During this time, researchers have been paying more attention to station-keeping by using electric propulsion. Oleson et al. investigated applying different electric propulsion systems for the north–south station-keeping [15]. Gomes and Prado used a hybrid optimal control approach to minimize the fuel consumption of the north–south station-keeping [16]. Gazzino et al. decomposed the minimum-fuel station-keeping problem into two steps: first, solve the control problem with an indirect method, and then deal with constraints [17]. They later presented further work in which a three-step decomposition method is employed to solve the same problem [18]. They also investigated solving the same problem via linear integer programming [19,20]. De Bruijn et al. linearized the relative element model in the geostationary station-keeping problem and came up with a method for determining station-keeping maneuvers by using convex optimization based on this model [21]. Min et al. discussed the influence of the configuration of electric thrusters on station-keeping and momentum-dumping [22]. Meanwhile, how to improve the autonomy of the on-board control system is also a focus. Guelmen creatively solved this problem by using a closed-loop control law, which requires the satellite’s real-time position and velocity information [23]. Because real-time GNSS determination is not always available for GEO satellites [24] and the output does not enforce constraints, deploying this method in the application would be challenging. Yang and Li studied the autonomous station-keeping control method for GEO satellites and developed the corresponding software [25]. Park et al. also proposed an autonomous station-keeping system based on the real-time determination [26].
In addition to the above studies, as a promising candidate for solving the station-keeping problem of electric-propulsion GEO satellites, the application of the model-predictive control (MPC) in this area has made significant progress in recent years. Being good at coordinating multiple inputs/outputs and handling constraints [27], the MPC is capable of optimizing the station-keeping performance while considering various constraints in the actual application. Walsh et al. studied applying the MPC for station-keeping and momentum management [28,29,30]. In their work, the controller used a prediction model based on the linearization of the orbital and attitude dynamics of GEO satellites to optimize the fuel consumption while enforcing the constraints caused by the thruster configuration. They chose the classical configuration in which four electric thrusters are fixed on the anti-nadir face of the satellite. Their result, given by the annual velocity increment, is Δ V AN = 78.67 m / s . In another work of their team, Weiss et al. deployed a six-thruster configuration to save fuel consumption [31] and achieved a better result of Δ V AN = 59.4 m / s . More thrusters means higher costs and lower reliability. Therefore, Caverly et al. proposed a new configuration where four thrusters are mounted on two assemblies, and they obtained Δ V AN = 66.8 m / s [32]. Although this result is higher than that in Ref. [28], it is still very competitive because they considered the on–off nature of the thruster. Unlike these three studies, Zou et al. considered the geostationary station-keeping problem with the MPC as an optimization problem based on the orbital elements instead of the position and velocity [33]. In their study, they combined the linearized relative element model proposed in Ref. [20] with the MPC and obtained Δ V AN = 72.3 m / s .
In this paper, in contrast to the previous studies, we apply a new configuration where thrusters are mounted on two symmetrical robotic arms. This configuration was first proposed by the Airbus company and had been applied in a space mission with the launch of EUTELSAT 172B in June 2017 [34]. They named it the Deployable Thruster Module Assembly (DTMA) and applied it to the Airbus Eurostar E3000 communication satellite platform. More information about this platform can be found in their introduction [35]. This platform includes a versatile system that enables the thruster to point in the optimal direction, benefiting both orbital transfer and station-keeping missions. Although our study is not based on the specific platform of the Airbus Eurostar E3000, we are inspired by the forward-looking idea of introducing robotic arms into the propulsion system. Therefore, in this paper, we discuss how to design a controller that can take advantage of the wider range of the torque-free thrust direction brought by the robotic arm. The detailed configuration applied in this paper will be stated in a later section.
In the previous MPC-related works, the researchers solved the station-keeping problem within a short prediction horizon based on linear dynamics. Then they quantized the optimal solution to enforce the constraints of the propulsion system, which could avoid handling the nonlinear constraints directly. In this way, the thrust direction will not be regarded as control variables in the control system. However, we believe that changing the thrust direction is critical for saving the fuel consumption of station-keeping. This will make the thrust direction subject to the nonlinear constraints caused by the pitch angle and the yaw angle of the thruster. Therefore, we have to deal with the nonlinear model predictive control (NLMPC) problem if introducing the MPC. The general approach to solving this kind of nonlinear optimization problem is to first find an initial value, and then use the optimization algorithm to obtain the solution, which is why we have to develop an open-loop solver first.
Therefore, in our study, we consider the approach to figuring out the station-keeping problem as three steps: first, we develop a state predictor (SP) to provide the prediction of the orbital states of GEO satellites; then based on our analysis of the GEO dynamics and the thruster configuration, we propose the quick feedforward controller (QFFC) as the open-loop solver, which could work out the solution for the current control cycle. Finally, we combine the SP and the QFFC with the NLMPC problem to form the fuel-optimal model predictive controller (FOMPC). The FOMPC applies a simulated annealing (SA) algorithm to find the optimal control solution with the initial value given by the QFFC. The FOMPC achieves less fuel consumption by optimizing the velocity increment of the north–south station-keeping, considering the station-keeping objective and the thruster configuration as constraints. Simulations show that both the QFFC and the FOMPC could achieve stable station-keeping. The QFFC has minor computational requirements and high accuracy for the station-keeping, while the FOMPC could achieve much less fuel consumption.
This paper is organized as follows. First, the electric thruster configuration and the dynamics are stated in Section 2. Next, in Section 3 and Section 4, we introduce the geostationary station-keeping problem and the proposed controllers in this paper. Then, we present several numerical simulations in Section 5 to evaluate the performance and robustness of the controllers. Finally, we draw conclusions in Section 6.

2. Thruster Configuration and Dynamics

In this paper, we assume a GEO communication satellite equipped with robotic arms whose nominal longitude and latitude of the sub-satellite point are λ T and 0 , and the main objective of station-keeping is to maintain the longitude and the latitude of the sub-satellite point. The thruster configuration and the nonlinear orbital dynamics we apply in this station-keeping problem will be stated in this section.

2.1. Thruster Configuration

As we have stated before, the thruster configuration in this paper is characterized by two symmetrical robotic arms placed at the anti-nadir face of the satellite. Thrusters are installed on the platform at the end of two robotic arms, providing the motorization to enable a wider range of thrust direction. Figure 1 shows the model of a GEO satellite that applies this configuration in the RSW (radial, along-track, cross-track) coordinate system, where C and C E are the centers of mass of the satellite and the Earth, respectively. and P represents the sub-satellite point on the ground, whose longitude and latitude are supposed to be λ T and 0 . The X-axis is in the orbital plane and points to the anti-nadir direction, while the Z-axis is perpendicular to the orbital plane. The robotic arms are installed on the two corners of the satellite’s anti-nadir face symmetrically. Due to the requirement of the antenna direction, the attitude of the satellite shown in Figure 1 is fixed during the working of the communication module.
With the help of the robotic arm, the thrusters are capable of flexible motion. The thrust direction in this paper is defined as the torque-free direction that passes through the center of mass of the satellite. Because momentum management is not involved in this paper, we do not consider the inside mechanics of the robotic arm. The thrust direction could be given by the yaw angle α and the pitch angle β , as shown in Figure 1. α represents the angle between the projection of the thrust direction in the X–Y plane and the +X-axis, while β represents the angle between the thrust direction and the X–Y plane. In this paper, we generally assume that thruster 1 could only work within ( + + + ) and ( + + ) space, and thruster 2 could only work within ( + + ) and ( + ) space. Thus, we could simplify the range of the feasible thrust direction to several inequality constraints, which could be given by
α min α 1 α max α max α 2 α min β min β 1 β max β max β 2 β min
where α 1 , β 1 , α 2 , β 2 are corresponding to thruster 1 and thruster 2, respectively, and α max , α min , β max , β min are nonnegative parameters. We will discuss the impact of these parameters in Section 5.

2.2. Dynamics

The acceleration of the GEO satellite a sat mainly consists of two parts, which could be written as
a sat = a thr + a p
where a thr is the acceleration generated by the thruster, and a p represents the perturbation acceleration.
GEO satellites are mainly affected by the solar and lunar gravitational attraction, the non-spherical perturbation of the Earth, and the solar radiation pressure [36]. Furthermore, the solar and lunar gravitational attraction not only changes the satellite’s inclination, but its components in the orbital plane also have a considerable impact on the semi-major axis and the eccentricity vector, making it the most important item among all perturbations. For the non-spherical perturbation of the Earth, we only consider the J 2 perturbation and ignore the remaining higher-order terms. The solar radiation pressure depends on many factors, such as the reflection coefficient and the surface-to-mass ratio of the satellite, and it does not have much influence on the orbital elements compared to the former two perturbations. Nevertheless, the solar radiation pressure is generally considered when studying the geostationary station-keeping problem due to its impact on the attitude of the satellite. Although this paper does not involve momentum management, we take the perturbation of solar radiation pressure into account to be consistent with other studies. Thus, these perturbation accelerations could be expressed as [31]
a p = a sun + a moon + a J 2 + a srp
where a sun , a moon , a J 2 , and a srp are the perturbations of the solar and lunar gravitational attraction, the J 2 perturbation, and the solar radiation pressure, respectively. Their analytic expressions are [37]
a q = μ q R q Δ q 3 r q r q 3 ( q represents the sun or the moon ) a J 2 = 3 μ J 2 r E 2 2 r 5 5 r · z ^ ECI r 2 1 r 2 r · z ^ ECI z ^ ECI a srp = ρ sun 1 + η S e m r sun 2 Δ sun 2 R ^ sun
where μ , μ sun , and μ moon are the gravitational constants of the Earth, the sun, and the moon, R q is the position vector of q relative to the satellite, r q is the position vector of q relative to the Earth, r is the position vector of the satellite, J 2 is a coefficient in the non-spherical perturbation model, z ^ ECI is the unit vector of the Z-axis in the Earth-centered inertial (ECI) coordinate system, ρ sun is the solar radiation pressure around the Earth, and η and S e m are the reflection coefficient and the effective surface-to-mass ratio of the satellite. a sat and all its components are in the ECI coordinate system, which would be more useful if transformed into the RSW coordinate system as a r , a t , a n T .
In order to avoid singularity, we adopt the modified equinoctial elements instead of the Kepler orbital elements in our study. The modified equinoctial elements are given by [38]
p = a 1 e 2 f = e cos ω + Ω g = e sin ω + Ω h = tan i / 2 cos Ω k = tan i / 2 sin Ω L = Ω + ω + θ
in which a , e , i , Ω , ω belong to the Kepler orbital elements, and θ is the true anomaly. The high-precision propagation of p , f , g , h , k , L is based on the Gauss variational equations [38], which could be written as
p ˙ = p μ 2 p w a t f ˙ = p μ a r sin L + 1 + w cos L + f a t w h sin L k cos L g · a n w g ˙ = p μ a r cos L + 1 + w sin L + g a t w + h sin L k cos L f · a n w h ˙ = p μ s 2 a n 2 w cos L k ˙ = p μ s 2 a n 2 w sin L L ˙ = μ p w p 2 + 1 w p μ h sin L k cos L a n
where w and s 2 are defined as
w = 1 + f cos L + g sin L s 2 = 1 + h 2 + k 2 .

3. Geostationary Station-Keeping Problem

In this section, we state the objective of geostationary station-keeping first. Then, combined with the thruster configuration mentioned in the previous section, we introduce some constraints into this problem. Finally, we present the control variables of the station-keeping problem and their relationship with the control effect.

3.1. Control Objectives

The primary objective of geostationary station-keeping is to prevent the longitude and the latitude of the GEO satellite from drifting too far. However, in addition to this primary objective, it is also necessary to control the eccentricity vector e x , e y T = [ f , g ] T , the inclination vector i x , i y T = [ h , k ] T , and the drift speed of longitude λ ˙ , which are important for the stability of station-keeping. Besides, some specific scenarios, such as the multi-satellite collocation, also have requirements for the eccentricity and inclination vectors [3]. Furthermore, due to the perturbations λ , e x , and e y have strong short-period oscillations, which makes it inconvenient for us to analyze their deviations. Therefore, we choose the mean elements to describe the orbital state of the GEO satellite, which could be written as
x = λ ¯ λ ¯ ˙ e ¯ x e ¯ y i ¯ x i ¯ y T .
The mean elements are obtained by the averaging operator [39]
ε ¯ = 1 T s 0 T s ε d t = 1 2 π 0 2 π ε d M
where T s is the orbital cycle of the GEO satellite, and ε ¯ donates the mean elements in Equation (8). Meanwhile, the objectives of the mean element are given by
x T = λ T λ ˙ T e x , T e y , T i x , T i y , T T .
Thus, the deviation δ x is
δ x = δ λ δ λ ˙ δ e x δ e y δ i x δ i y T = λ T λ ¯ λ ˙ T λ ¯ ˙ e x , T e ¯ x e y , T e ¯ y i x , T i ¯ x i y , T i ¯ y T .
The primary objective of station-keeping is to compensate for δ x so as to ensure that x = x T at the end of each control cycle, thus giving the equality constraints of the station-keeping objective.

3.2. Constraints

In this paper, we consider many constraints while solving the optimization problem instead of enforcing them after obtaining the optimal solution. Some of these constraints are imposed on the satellite, while others are imposed on the station-keeping process. The satellite model equipped with robotic arms that we study in our work is a very specific model. Therefore, we expect that the constraints of the problem could also be close to those in the real application. These constraints could be summarized as follows.
  • The thrust direction is limited by the robotic arms, as stated in Section 2.1.
  • Although the thrust direction could be adjusted with the robotic arm, the adjustment is supposed to be completed before the maneuver. The robotic arms are locked during the working of thrusters.
  • Because station-keeping is carried out at a fixed cycle, the length of each control cycle determines how far the satellite will drift during this time, making it essential to improve the control accuracy. Meanwhile, the control cycle should be the integer multiples of the orbital cycle to simplify the control problem. Therefore, the control cycle is assumed to be one orbital cycle in our work.
  • Only one thruster can work at the same time due to the on-board power limitation.
  • The reliable ignition times of the thruster are limited. In order to reduce the ignition times, each thruster only ignites once in one control cycle, and each ignition corresponds to one maneuver. For one control cycle, maneuver 1 is performed by thruster 1, and maneuver 2 is performed by thruster 2.
  • The only exception to constraint 5 is that the satellite happens to enter the Earth’s shadow while the thruster is still working. Due to the lack of energy, the working of the thruster will be suspended until the satellite leaves the shadow region. In this case, one thruster actually carries out two maneuvers before and after entering the shadow region, but they will be regarded as one maneuver in the optimization.
  • To reduce the ignition times, the minimum working time for each maneuver is 30 s.

3.3. Control Solution

Based on the above description, a complete control solution S for one control cycle should involve eight variables corresponding to the two maneuvers. S could be written as
S = [ α 1 , β 1 , α 2 , β 2 , l 1 , l 2 , V 1 , V 2 ] T
where α 1 , α 2 , β 1 , and β 2 are the yaw angles and pitch angles of two thrusters, l 1 and l 2 are the satellite’s right ascensions when conducting the two maneuvers, and V 1 and V 2 are the velocity increments of the two maneuvers. Because the maneuver usually does not last too long, l 1 and l 2 could be regarded as the value at the middle, and V 1 , V 2 could be approximate to [40]
V i = F T m D i ( i = 1 , 2 )
in which F T is the thrust magnitude generated by the thruster, D i is the length of the i-th maneuver, and m is the mass of the satellite at the beginning of the control cycle, neglecting small mass changes during the maneuver. The relationship between S and the velocity increments v 1 , v 2 in the RSW coordinate system could be written as
v i = v r i v t i v n i = V i cos β i cos α i cos β i sin α i sin β i ( i = 1 , 2 ) .
The changes of the mean elements caused by the two maneuvers are defined as Δ x 1 and Δ x 2 , which are given by
Δ x i = Δ λ i Δ λ ˙ i Δ e x i Δ e y i Δ i x i Δ i y i T = B i v i ( i = 1 , 2 )
in which B i is derived from the approximate form of the Gaussian equations [39]
a ˙ = 2 R s V s a t
λ ˙ = n Ω E 2 V s a r
e ˙ x = 1 V s a r sin l + 2 a t cos l
e ˙ y = 1 V s a r cos l + 2 a t sin l
i ˙ x = cos l V s a n
i ˙ y = sin l V s a n
where R s and V s are the nominal radius and velocity of GEO satellites, n = μ a 3 and l = Ω + ω + θ are the mean angular velocity and the right ascension, and Ω E is the angular rotation rate of the Earth. It could be derived from Equation (16b) that
Δ λ i = 2 V s v r i + Δ n i · T i
Δ λ ˙ i = Δ n i
where Δ n i is the change of n caused by the i-th maneuver, which is given by
Δ n i = d n d a · Δ a i = 3 2 μ a 5 · 2 R s V s v t 3 R s v t
and T i is the remaining time of the current control cycle after the i-th maneuver. T i could be approximate to
T i c i T s .
The parameter c i is given by
c i = 1 l i l start 2 π , l i l start l start l i 2 π , l i < l start , 0 l i , l start < 2 π
in which l start is the satellite’s right ascension when the current control cycle starts. In conclusion, B i could be written as
B i = 2 V s 3 c i V s 0 0 3 R s 0 sin l i V s 2 cos l i V s 0 cos l i V s 2 sin l i V s 0 0 0 cos l i 2 V s 0 0 sin l i 2 V s = .
Hence, defining the changes of mean elements caused by maneuvers in one control cycle as the control effect Δ x s , we could establish the relationship between S and Δ x s as
Δ x s = Δ x 1 + Δ x 2 = B 1 v 1 + B 2 v 2 = F S .

4. Station-Keeping Controller Design

In this section, we present the controller design for the geostationary station-keeping problem in this paper, which consists of three parts: the state predictor (SP), the quick feedforward controller (QFFC), and the fuel-optimal model predictive controller (FOMPC).

4.1. State Predictor

Let us define the mean elements at the beginning of the first control cycle as the initial state x 0 and the prediction at the end of the k-th control cycle as x k . This can be written as
x k = λ ¯ k λ ¯ ˙ k e ¯ x , k e ¯ y , k i ¯ x , k i ¯ y , k T ( k = 0 , 1 , 2 , . . . ) .
The primary function of the SP is to obtain the prediction x k , given the initial state x 0 and the control solutions of former control cycles S 1 , S 2 , . . . , S k .
Generally, the element changes caused by station-keeping maneuvers and natural drifts are small within several control cycles, making them independent of each other. Thus, the x k can be obtained by
x k = x k 1 + Δ x s , k + Δ x d , k
in which Δ x s , k and Δ x d , k are the element changes caused by maneuvers and natural drifts, respectively. Δ x s , k can be obtained through Equation (23), while Δ x d , k can be written as
Δ x d , k = Δ λ d , k Δ λ ˙ d , k Δ e x d , k Δ e y d , k Δ i x d , k Δ i y d , k T .
Δ x d , k could be easily obtained in advance through the high-precision orbit propagation. The high-precision orbit propagation is completed in the preparation of the SP and only needs to be conducted once, causing a minor computational burden. Besides, Δ λ d , k could be approximate to
Δ λ d , k λ ¯ ˙ k 1 T s .

4.2. Quick Feedforward Controller

With the help of the SP, we present the quick feedforward controller (QFFC). The primary function of the QFFC is to quickly work out the control solution for the current control cycle based on the prediction given by the SP. Because the control objective is to compensate for the deviation of the element, it is required that x k = x T . Therefore, according to Equation (25), the effect of the control solution is supposed to fulfill
Δ x s , k = x T x k 1 Δ x d , k .
Defining δ x k = x T x k 1 Δ x d , k as the predictive deviation of the k-th control cycle, we can obtain the basic rule of the feedforward control as
δ x k = Δ x s , k = B 1 , k v 1 , k + B 2 , k v 2 , k = F S k
in which B 1 , k , B 2 , k , v 1 , k , and v 2 , k are given by Equations (14) and (22).
For S k = [ α 1 , k , β 1 , k , α 2 , k , β 2 , k , l 1 , k , l 2 , k , V 1 , k , V 2 , k ] T , we make an important simplification first. The last two rows of Equation (29) can be written as
δ i x , k δ i y , k = v n 1 , k 2 V s cos l 1 k sin l 1 k + v n 2 , k 2 V s cos l 2 , k sin l 2 , k .
Because the north–south station-keeping dominates the required velocity increment, it must be carried out as efficiently as possible to save fuel consumption. According to Equation (30), [ cos l 1 , k , sin l 1 , k ] T and [ cos l 2 , k , sin l 2 , k ] T have to be in the same or opposite direction as [ δ i x , k , δ i y , k ] to maximize the efficiency of north–south station-keeping, meaning that either l 2 , k = l 1 , k or l 2 , k = l 1 , k + π . If l 2 , k = l 1 , k , the two maneuvers will be equivalent to one larger maneuver, contrary to our original intention of setting two maneuvers. Therefore, based on l 2 , k = l 1 , k + π , we could obtain that
l 2 , k = l 1 , k + π = arctan δ i y , k δ i x , k .
This simplification would reduce the number of variables in the control solution from eight to six.
Then Equation (29) could be rewritten as
V L , k = v r 1 , k + v r 2 , k + 3 2 ( c 1 , k v t 1 , k + c 2 , k v t 2 , k )
V DT , k = v t 1 , k + v t 2 , k
V R , k = v r 1 , k v r 2 , k
V ET , k = v t 1 , k v t 2 , k
V N , k = v n 1 , k + v n 2 , k
among which V L , k , V DT , k , V R , k , V ET , k , V N , k are defined as
V L , k = V s 2 δ λ k
V DT , k = R s 3 δ λ ˙ k
V R , k = V s δ e x , k sin l 1 , k δ e y , k cos l 1 , k
V ET , k = V s 2 δ e x , k cos l 1 , k + δ e y , k sin l 1 , k
V N , k = 2 V s cos l 1 δ i x , k = 2 V s sin l 1 δ i y , k .
So far, if the control system is not subject to any constraint, it would be easy to obtain the solution by solving Equations (32a)–(32e). It is worth noting that there is usually more than one solution for Equations (32a)–(32e), among which we choose the one that minimizes V 1 , k + V 2 , k as the control solution.
Furthermore, because the thruster cannot face the nadir direction as stated in Section 2.1, it is obvious that
v r 1 , k 0 , v r 2 , k 0 .
Then by combining Equation (34) with Equations (32a)–(32d), we could obtain
V L , k 3 4 c 1 , k + c 2 , k V DT , k 3 4 c 1 , k c 2 , k V ET , k + V R , k 0 .
Equation (35) is a necessary condition for Equations (32a)–(32e) to have a solution. If the V L , k , V DT , k , V R , k , V ET , k are all determined, this condition may not be matched. To solve this problem, we have to change some of them to make sure the condition of Equation (35) is satisfied. Changing V L , k , V R , k , V ET , k means changing δ λ k , δ e x , k , δ e y , k , which will directly affect the stability of station-keeping, while changing V DT , k (i.e., changing δ λ ˙ k ) will not. Therefore, increasing V DT , k until we obtain a solution is the only option under these circumstances.
Although V L , k has been determined for the current control cycle, we can still change V L , k + 1 indirectly by changing V DT , k , making it easier for the condition of Equation (35) to be matched in the next control cycle. The condition of Equation (35) in the next control cycle can be written as
V L , k + 1 3 4 c 1 , k + 1 + c 2 , k + 1 V DT , k + 1 + 3 4 c 1 , k + 1 c 2 , k + 1 V ET , k + 1 V R , k + 1
in which V ET , k + 1 and V R , k + 1 are determined by the prediction of δ x k + 1 , and V DT , k + 1 is related to later control cycles. Thus, we assume that V DT , k + 1 = 0 and obtain the estimate of the required V L , k + 1 in the next control cycle as
V ˜ L , k + 1 = 3 4 c 1 , k + 1 c 2 , k + 1 V ET , k + 1 V R , k + 1 .
Because
λ ¯ ˙ k Δ λ d , k + 1 T s = λ T λ ¯ k δ λ k + 1 T s = δ λ k + 1 T s = 2 R s V ˜ L , k + 1 ,
the estimate of the required V DT , k can be given by
V ˜ DT , k = R s 3 δ λ ˙ k = R s 3 λ ¯ ˙ k λ ¯ ˙ k 1 Δ λ ˙ d , k = 2 3 V ˜ L , k + 1 + R s 3 λ ¯ ˙ k 1 + Δ λ ˙ d , k .
The intention of this step is not to guarantee the satisfaction of the condition of Equation (35). In fact, it would not help the situation of the current control cycle. If V ˜ DT , k cannot satisfy the condition of Equation (35), we still need to increase V DT , k based on V ˜ DT , k until we obtain a solution. Changing V L , k + 1 indirectly could avoid a large adjustment to V DT , k + 1 in the next control cycle. On the whole, this will smooth the variation of λ ¯ ˙ , thereby enhancing the stability of the station-keeping and slightly reducing fuel consumption.
In summary, by establishing the relationship between the control solution and the control objective, we design the QFFC based on the feedforward control. The QFFC mainly has two roles: an open-loop solver in the following FOMPC or an autonomous on-board controller. We will evaluate the performance of the QFFC in Section 5.

4.3. Fuel-Optimal Model Predictive Controller

The QFFC stated above obtains the solution by adjusting V DT . Not optimal though the V DT is, it has little influence on the overall fuel consumption due to its small proportion of the required velocity increment. Actually, the required velocity increment is dominated by north–south station-keeping. Therefore, to effectively optimize the fuel consumption, the controller needs to adjust the V N of each control cycle according to the prediction given by the SP. On the other hand, the analysis of the condition of Equation (35) is only concerned with the in-plane components of the velocity increment. Under extreme circumstances, the QFFC may not be able to give the solution, which has to be fixed by adjusting V N . Adjusting V N means we just need to ensure that the deviation of the inclination vector does not exceed the threshold Δ i max at the end of each control cycle instead of requiring i ¯ x , i ¯ y T = i x , T , i y , T T strictly. Defining the value of V N applied in the k-th control cycle as u k , we could regard the process to obtain the control solution by the QFFC as a function of u k , which can be written as
S k = Q u k δ x = δ x k , if solution exists 0 , if solution does not exist .
Therefore, the nonlinear model predictive control (NLMPC) problem can be formulated as
min U h 1 2 k = 1 N ( V 1 , k | h 2 + V 2 , k | h 2 ) s . t . x k | h = x k 1 | h + Δ x s , k | h + Δ x d , k | h , x 0 | h = x h , λ ¯ k | h = λ T , e ¯ x , k | h = e x , T , e ¯ y , k | h = e y , T , i x , T i ¯ x , k | h 2 + i y , T i ¯ y , k | h 2 Δ i max , α min α 1 , k | h α max , α max α 2 , k | h α min , β min β 1 , k | h β max , β max β 2 , k | h β min , V 1 , k | h 0 , V 2 , k | h 0
where h is the number of the current control cycle, N is the prediction horizon, ς k | h denotes the prediction of ς for k control cycles ahead of h, U h = u 1 | h , u 2 | h , . . . , u N | h T is the sequence of u k | h within the prediction horizon, S k | h = [ α 1 , k | h , β 1 , k | h , α 2 , k | h , β 2 , k | h , V 1 , k | h , V 2 , k | h ] T are given by Equation (40). The optimal solution of Equation (41) could be written as U h * , among which we only apply u 1 | h * to the current control cycle h following the rule of the MPC.
The prediction horizon N determines the number of control cycles involved in a single optimization. When N = 1 , Equation (41) becomes a greedy algorithm for minimizing the fuel consumption of the current cycle; when N = , the solution of Equation (41) is the actual global optimum of the station-keeping problem. Therefore, increasing N benefits the performance of the controller in theory. Whereas with the growth of the prediction horizon, the SP cannot keep the prediction accuracy, and the computational requirement increases fast as well. Therefore, the optimal value of N exists. We will verify this in the next section.
There are several practical algorithms to solve the nonlinear optimization problem [41], among which we employ a simulated annealing (SA) algorithm in the optimization. The algorithm was first proposed to solve the combinational optimization problem by simulating the solid annealing phenomenon [42]. The SA has good robustness and parallelism, but it also turns out to be sensitive to the initial value. In our case, the SA is used to find the optimal solution sequence U h * . Based on Equation (41), we define the evaluation function as
J U h   = 1 2 k = 1 N V 1 , k | h 2 + V 2 , k | h 2 + ξ k i x , T i ¯ x , k | h 2 + i y , T i ¯ y , k | h 2 , ξ k = 0 , i x , T i ¯ x , k | h 2 + i y , T i ¯ y , k | h 2 Δ i max P , i x , T i ¯ x , k | h 2 + i y , T i ¯ y , k | h 2 > Δ i max
where ξ k and P are penalty parameters to ensure that U h , which fails to enforce the constraint of the inclination vector, has poor performance. The initial value U h 0 of the optimization is given by the QFFC. The algorithm explores the solution space by generating disturbances. In the end, as the temperature continues to drop, and the result of the SA converges to the optimal solution U h * . Applying u 1 | h * to the current control cycle h, we obtain the fuel-optimal control solution. The performance of the FOMPC will be evaluated in Section 5.

5. Numerical Simulations

In this section, we present the results of several numerical simulations to test the QFFC and the FOMPC. Then, by setting different constraints and introducing the attitude and determination errors, we assess the robustness of the two controllers. To carry out simulations, we assume a testing environment with a GEO satellite of initial mass m 0 = 5000 kg . The satellite is equipped with two identical thrusters, which could generate the constant thrust of F T = 500 mN and have the specific impulse of I sp = 2250 s . We will analyze the influence of different constraints and parameters in some simulations, while in other simulations, they are set to their nominal values, which are β max = 70 , β min = 0 , α max = 60 , α min = 30 , Δ i max = 8 × 10 5 , and N = 6 . The dynamics and perturbations are stated in Section 2.2. Let the coefficients of the solar radiation pressure in Equation (4) be ρ sun = 4.5605 × 10 6 N / m 2 , η = 0.6 , and S = 200 m 2 . The controllers are based on the mean elements, whereas the simulations are conducted by the high-precision orbit propagator with the osculating modified equinoctial elements. The start time of the mission is set to 12:00 on 1 January 2020, and the initial osculating elements of the GEO satellite are
p , f , g , h , k , L = 42167 , 0 , 1 × 10 4 , 0 , 5 × 10 5 , 0.289893 .
The objective elements are given by
λ T , e x , T , e y , T , i x , T , i y , T = 96 , 0 , 1 × 10 4 , 0 , 3 × 10 5 .
The algorithms are programmed in Fortran 90, and the testing platform is Intel® Core™ i7-9700K 3.6 GHz. To be consistent with the previous studies, we measure the fuel consumption of station-keeping with the annual velocity increment V ¯ AN .

5.1. Results

First, we carry out the one-year simulation for the two controllers. Figure 2 and Figure 3 present the time history of the osculating elements in the two simulations. Figure 4, Figure 5 and Figure 6 compare the trajectories in the ( δ λ , δ ϕ ) , ( δ e x , δ e y ) , and ( δ i x , δ i y ) plane, in which the red dashed line represents the trajectory in the uncontrolled situation. It’s obvious that without station-keeping, the satellite will quickly drift away from the objective. The results of the two simulations show that both the QFFC and the FOMPC could achieve long-term stable station-keeping with the V ¯ AN of 60.2 m / s / year and 54.2 m / s / year , respectively. The number of maneuvers during 365 orbits of the QFFC and the FOMPC are 728 and 721, suggesting that each thruster performs approximately one maneuver per orbit. The QFFC and the FOMPC have almost the same performance in controlling the longitude and the eccentricity vector. According to Figure 4 and Figure 5, the control accuracy for the longitude and the eccentricity vector is about ± 0.013 and 1.1 × 10 4 , respectively. The main difference between the two simulations is that the QFFC shows higher accuracy for controlling the inclination vector. That is because the QFFC requires the inclination vector to be fully adjusted to the objective at the end of each control cycle, while the FOMPC allows a certain deviation. Besides, Figure 3 shows that the inclination vector presents a more apparent periodic variation under the control of the FOMPC. To figure out the reason, we give Figure 7, which shows the half-year time history of the north–south velocity increment Δ V NS and Figure 8, which shows that of the average pitch angle β ¯ . Δ V NS and β ¯ are given by
Δ V NS = V 1 sin β 1 + V 2 sin β 2
β ¯ = arcsin ( Δ V NS V 1 + V 2 ) .
Compared with the blue curve representing the FOMPC in Figure 7, the red one representing the QFFC shows a larger variation, reflecting the impact of the solar and lunar gravitational attraction. Meanwhile, according to Figure 8, β ¯ reaches its ceiling of 70 most of the time in both simulations, indicating that the east–west velocity increment is usually forced to increase to adapt to the north–south velocity increment due to the constraint of the pitch angle. Therefore, by reducing the variation of Δ V NS , the above situation may be avoided, thereby improving the efficiency of station-keeping, explaining how the FOMPC manages to save fuel consumption. However, this will allow the inclination vector to vary with the solar and lunar gravitational attraction, reducing the control accuracy of the inclination vector, as shown in Figure 6. The average CPU times per orbit of the QFFC and the FOMPC are 0.09 s and 57.4 s, respectively. Although the FOMPC has higher computational requirements, considering further code optimization and the upgrade of the computing system on the satellite, the on-board application of the FOMPC is possible.
Next, we conduct the following simulations concerning the impact of different parameters and constraints. Due to the characteristic of the optimization algorithm, the results of each simulation will be slightly different, even if the parameters are the same. Thus, we apply the Monte Carlo simulation method [43], in which we repeat the same simulation 20 times and take the average annual velocity increment as the V ¯ AN . Figure 9, Figure 10, Figure 11 and Figure 12 present the simulations for the constraints of the yaw angle α and the pitch angle β . When we change one of these parameters, the others are assumed to be β max = 90 , β min = 0 , α max = 90 , α min = 90 . The results show that the effects of β max , β min , α max , α min on the QFFC and the FOMPC are almost the same. β max has the largest impact on the efficiency of station-keeping. As β max increases from 45 to 90 , the V ¯ AN drops from 78.2 m / s / year to 57.5 m / s / year for the QFFC and 73.7 m / s / year to 51.4 m / s / year for the FOMPC. The drop is more apparent when β max is less than 80 , indicating that β max will remain a major limitation on the station-keeping efficiency until it reaches 80 . On the contrary, according to Figure 10, β min has almost no influence on the V ¯ AN until it reaches 40 . According to Figure 11 and Figure 12, α max and α min have the similar impact due to the symmetry. The restricted areas in the two figures represent the values of α max and α min with which the controllers could not operate stably. Outside the areas, the influences of α max and α min declines fast and becomes negligible after α max and α min reach 45 .
For the parameters N and Δ i max in the FOMPC, Figure 13 proves the existence of the optimal value of N, which are 3, 6, and 7 under weak, normal, and strong constraints. The detailed information of these three levels can be found in Table 1. It is obvious that the optimal value of N goes up as the constraints get stronger, for the reason that strong constraints make the optimization more necessary. Figure 14 indicates that Δ i max significantly impacts the V ¯ AN because it determines the degree of optimization.
Finally we demonstrate the robustness of the two controllers by taking the attitude and determination errors into account. We assume that these disturbances follow the Gaussian distribution as
χ N χ true , σ χ 2
where χ donates the yaw angle α , the pitch angle β , the position vector r , and the velocity vector v , respectively, χ real is the true value, and σ χ determines the error distribution range (EDR). We assume the nominal error distribution range (NEDR) of the attitude and determination as
3 σ α , β = 0 . 05 3 σ r = 0.3 m , 3 σ v = 5 × 10 5 m / s ,
which represents a high-accuracy level considering the present technology. Magnifying NEDR by H χ times and applying them in the algorithm, we perform similar Monte Carlo simulations to investigate the influence of the errors of different types and different magnitudes. As shown in Figure 15, the attitude error has almost no effect on the two controllers. However, Figure 16 suggests that raising the EDR of the determination results in a significant increase of the V ¯ AN . This is because the determination error directly affects the accuracy of the predictive model. For both the QFFC and the FOMPC, the station-keeping becomes unstable after H r , v exceeds 800, indicating that the maximum EDR of the determination is 3 σ r = 240 m , 3 σ v = 0.04 m / s , which is closed to the accuracy of the on-board real-time GNSS determination [44]. In Figure 17, we also present the trajectories in the ( δ λ , δ ϕ ) plane with the maximum EDR of the determination. Compared with Figure 4, the control accuracy of the longitude has risen to ± 0.025 for the QFFC and ± 0.030 for the FOMPC. However, that of the latitude hardly changes, indicating that the stability of station-keeping for both the QFFC and the FOMPC is guaranteed under the condition of low determination accuracy.

5.2. Discussion

In the last part of this paper, we evaluate the performance of the QFFC and the FOMPC, considering the parameters and constraints of different values. Simulation results verify the validity of the QFFC and the FOMPC. Among all these parameters, β max has the most considerable influence on fuel consumption; α max , α min , and β min affect fuel consumption significantly only if the constraints are too strong. For the FOMPC, Δ i max determines the degree of optimization; the predictive horizon N has a certain impact, and the optimal value of N depends on how strong the constraints are. Robustness verification proves that both controllers could handle the attitude and determination errors of different magnitudes, while the rise of the determination error will cause extra fuel consumption. Here, we present Table 2, in which we summarize our results under normal constraints and compare them with those in the previous studies. The results prove that our approaches are effective.

6. Conclusions

Powered by advanced electric propulsion technology, the geostationary station-keeping problem using electric propulsion has attracted much attention. In this paper, we apply a novel thruster configuration involving two robotic arms. The QFFC and the FOMPC are proposed, and an NLMPC problem is solved. Both the QFFC and the FOMPC are proven to be effective in controlling multiple elements of the GEO satellite. The QFFC is capable of carrying out high-accuracy station-keeping with minor computational requirements, whereas the FOMPC could significantly reduce fuel consumption. Further investigations suggest that the parameters and constraints affect fuel consumption differently, among which the ceiling of the pitch angle of the thrust direction is the most influential parameter. The robustness assessments of the two controllers are also presented, considering the attitude and determination errors of different magnitudes. It proves that the two controllers are robust, but the rise of the determination error will cause extra fuel consumption.
Our study suggests that introducing robotic arms could benefit geostationary station-keeping significantly because it enables a wider range of the thrust direction, but its impact on the momentum management of the satellite is not involved in this paper. Next, the control problem that combines station-keeping and momentum management based on this thruster configuration will be our focus. Although the real application scenarios may be different from our assumptions, we expect our study to provide references for the control system design of GEO missions in the future.

Author Contributions

Conceptualization, W.Z. and Q.P.; Writing—original draft, C.L.; Writing—review and editing, B.X. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Manned Space Project Foundation of China (ZYZS12-QW04) and the Basic Research Project (JCKY2020110C096).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lemmens, M. Geo-Information: Technologies, Applications and the Environment; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2011; Volume 5. [Google Scholar]
  2. Finch, M.J. Limited space: Allocating the geostationary orbit. Nw. J. Int’l L. Bus. 1985, 7, 788. [Google Scholar]
  3. Lee, B.S.; Lee, J.S.; Choi, K.H. Analysis of a station-keeping maneuver strategy for collocation of three geostationary satellites. Control Eng. Pract. 1999, 7, 1153–1161. [Google Scholar] [CrossRef]
  4. Shrivastava, S.K. Orbital perturbations and stationkeeping of communication satellites. J. Spacecr. Rocket. 1978, 15, 67–78. [Google Scholar] [CrossRef]
  5. Slavinskas, D.; Dabbaghi, H.; Benden, W.; Johnson, G. Efficient inclination control for geostationary satellites. J. Guid. Control Dyn. 1988, 11, 584–589. [Google Scholar] [CrossRef]
  6. Jahn, R.G. Electric propulsion. Am. Sci. 1964, 52, 207–217. [Google Scholar]
  7. Yang, G. Earth-moon trajectory optimization using solar electric propulsion. Chin. J. Aeronaut. 2007, 20, 452–463. [Google Scholar] [CrossRef] [Green Version]
  8. Yang, D.; Xu, B.; Zhang, L. Optimal low-thrust spiral trajectories using Lyapunov-based guidance. Acta Astronaut. 2016, 126, 275–285. [Google Scholar] [CrossRef]
  9. Takei, Y.; Saiki, T.; Yamamoto, Y.; Mimasu, Y.; Takeuchi, H.; Ikeda, H.; Ogawa, N.; Terui, F.; Ono, G.; Yoshikawa, K.; et al. Hayabusa2’s station-keeping operation in the proximity of the asteroid Ryugu. Astrodynamics 2020, 4, 349–375. [Google Scholar] [CrossRef]
  10. Feuerborn, S.A.; Perkins, J.; Neary, D.A. Finding a way: Boeing’s all electric propulsion satellite. In Proceedings of the 49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, San Jose, CA, USA, 14–17 July 2013; p. 4126. [Google Scholar]
  11. Jahn, R.G. Physics of Electric Propulsion; Courier Corporation: Honolulu, HI, USA, 2006. [Google Scholar]
  12. Hofer, R. High-specific impulse operation of the BPT-4000 Hall thruster for NASA science missions. In Proceedings of the 46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Nashville, TN, USA, 25–28 July 2010; p. 6623. [Google Scholar]
  13. Sankovic, J.M.; Hamley, J.A.; Haag, T.W. Performance evaluation of the Russian SPT-100 thruster at NASA LeRC. In Proceedings of the IEPC Conference, Seattle, WA, USA, 13–14 September 1993. [Google Scholar]
  14. Valentian, D.; Maslennikov, N. The PPS 1350 program. CP IEPC 1997, 97, 855–861. [Google Scholar]
  15. Oleson, S.R.; Myers, R.M.; Kluever, C.A.; Riehl, J.P.; Curran, F.M. Advanced propulsion for geostationary orbit insertion and north-south station keeping. J. Spacecr. Rocket. 1997, 34, 22–28. [Google Scholar] [CrossRef] [Green Version]
  16. Gomes, V.M.; Prado, A.F. Low-thrust out-of-plane orbital station-keeping maneuvers for satellites. Math. Probl. Eng. 2012, 2012, 532708. [Google Scholar] [CrossRef] [Green Version]
  17. Gazzino, C.; Arzelier, D.; Losa, D.; Louembet, C.; Pittet, C.; Cerri, L. Optimal control for minimum-fuel geostationary station keeping of satellites equipped with electric propulsion. IFAC-PapersOnLine 2016, 49, 379–384. [Google Scholar] [CrossRef]
  18. Gazzino, C.; Arzelier, D.; Cerri, L.; Losa, D.; Louembet, C.; Pittet, C. A three-step decomposition method for solving the minimum-fuel geostationary station keeping of satellites equipped with electric propulsion. Acta Astronaut. 2019, 158, 12–22. [Google Scholar] [CrossRef] [Green Version]
  19. Gazzino, C.; Louembet, C.; Arzelier, D.; Jozefowiez, N.; Losa, D.; Pittet, C.; Cerri, L. Integer programming for optimal control of geostationary station keeping of low-thrust satellites. IFAC-PapersOnLine 2017, 50, 8169–8174. [Google Scholar] [CrossRef]
  20. Gazzino, C.; Arzelier, D.; Louembet, C.; Cerri, L.; Pittet, C.; Losa, D. Long-term electric-propulsion geostationary station-keeping via integer programming. J. Guid. Control Dyn. 2019, 42, 976–991. [Google Scholar] [CrossRef]
  21. De Bruijn, F.J.; Theil, S.; Choukroun, D.; Gill, E. Geostationary satellite station-keeping using convex optimization. J. Guid. Control Dyn. 2016, 39, 605–616. [Google Scholar] [CrossRef] [Green Version]
  22. Min, W.; Qiang, L.; Xingang, L.; Ran, A. Electric Thrusters Configuration Strategy Study on Station-Keeping and Momentum Dumping of GEO Satellite. In Proceedings of the 2019 5th International Conference on Control Science and Systems Engineering (ICCSSE), Shanghai, China, 14–16 August 2019; pp. 115–121. [Google Scholar]
  23. Guelman, M.M. Geostationary satellites autonomous closed loop station keeping. Acta Astronaut. 2014, 97, 9–15. [Google Scholar] [CrossRef]
  24. Ruiz, J.L.; Frey, C.H. Geosynchronous satellite use of GPS. In Proceedings of the 18th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2005), Long Beach, CA, USA, 13–16 September 2005; pp. 1227–1232. [Google Scholar]
  25. Yang, W.; Li, S. A station-keeping control method for GEO spacecraft based on autonomous control architecture. Aerosp. Sci. Technol. 2015, 45, 462–475. [Google Scholar] [CrossRef]
  26. Park, B.K.; Tahk, M.J.; Bang, H.C.; Park, C.S.; Jin, J.H. A new approach to on-board stationkeeping of geo-satellites. Aerosp. Sci. Technol. 2005, 9, 722–731. [Google Scholar] [CrossRef]
  27. Garcia, C.E.; Prett, D.M.; Morari, M. Model predictive control: Theory and practice-A survey. Automatica 1989, 25, 335–348. [Google Scholar] [CrossRef]
  28. Weiss, A.; Di Cairano, S. Opportunities and potential of model predictive control for low-thrust spacecraft station-keeping and momentum-management. In Proceedings of the 2015 European Control Conference (ECC), Linz, Austria, 15–17 July 2015; pp. 1370–1375. [Google Scholar]
  29. Weiss, A.; Kalabić, U.; Di Cairano, S. Model predictive control for simultaneous station keeping and momentum management of low-thrust satellites. In Proceedings of the 2015 American Control Conference (ACC), Chicago, IL, USA, 1–3 July 2015; pp. 2305–2310. [Google Scholar]
  30. Walsh, A.; Di Cairano, S.; Weiss, A. MPC for coupled station keeping, attitude control, and momentum management of low-thrust geostationary satellites. In Proceedings of the 2016 American Control Conference (ACC), Boston, MA, USA, 6–8 July 2016; pp. 7408–7413. [Google Scholar]
  31. Weiss, A.; Kalabić, U.V.; Di Cairano, S. Station keeping and momentum management of low-thrust satellites using MPC. Aerosp. Sci. Technol. 2018, 76, 229–241. [Google Scholar] [CrossRef]
  32. Caverly, R.J.; Di Cairano, S.; Weiss, A. Electric Satellite Station Keeping, Attitude Control, and Momentum Management by MPC. IEEE Trans. Control Syst. Technol. 2020, 29, 1475–1489. [Google Scholar] [CrossRef]
  33. Zou, H.; Song, J.; Wang, J.; Zhang, L.; Huang, Y. Geostationary Station Keeping Using Relative Orbital Elements with Model Predictive Control. In Proceedings of the 2020 Chinese Control And Decision Conference (CCDC), Hefei, China, 22–24 August 2020; pp. 4437–4442. [Google Scholar]
  34. Boniface, C.; Charbonnier, J.M.; Lefebvre, L.; Leroi, V.; Lienart, T. An overview of electric propulsion activities at CNES. In Proceedings of the 35th International Electric Propulsion Conference, Atlanta, GA, USA, 8–12 October 2017. [Google Scholar]
  35. Sembély, X.; Wartelski, M.; Doubrère, P.; Deltour, B.; Cau, P.; Rochard, F. Design and Development of an Electric Propulsion Deployable Arm for Airbus Eurostar E3000 ComSat Platform. In Proceedings of the 35th International Electric Propulsion Conference, Atlanta, GA, USA, 8–12 October 2017. [Google Scholar]
  36. Curtis, H. Orbital Mechanics for Engineering Students; Butterworth-Heinemann: Oxford, UK, 2013. [Google Scholar]
  37. De Ruiter, A.H.; Damaren, C.; Forbes, J.R. Spacecraft Dynamics and Control: An Introduction; John Wiley & Sons: Hoboken, NJ, USA, 2012. [Google Scholar]
  38. Walker, M.; Ireland, B.; Owens, J. A set modified equinoctial orbit elements. Celest. Mech. 1985, 36, 409–419. [Google Scholar] [CrossRef]
  39. Li, L.; Zhang, J.; Zhao, S.; Qi, R.; Li, Y. Autonomous onboard estimation of mean orbital elements for geostationary electric-propulsion satellites. Aerosp. Sci. Technol. 2019, 94, 105369. [Google Scholar] [CrossRef]
  40. Boyd, I.D. Numerical modeling of spacecraft electric propulsion thrusters. Prog. Aerosp. Sci. 2005, 41, 669–687. [Google Scholar] [CrossRef]
  41. Bazaraa, M.S.; Sherali, H.D.; Shetty, C.M. Nonlinear Programming: Theory and Algorithms; John Wiley & Sons: Hoboken, NJ, USA, 2013. [Google Scholar]
  42. Sahin, K.H.; Ciric, A.R. A dual temperature simulated annealing approach for solving bilevel programming problems. Comput. Chem. Eng. 1998, 23, 11–25. [Google Scholar] [CrossRef]
  43. Mooney, C.Z. Monte Carlo Simulation; Sage: New York, NY, USA, 1997; Number 116. [Google Scholar]
  44. Xiong, Z.; Qiao, L.; Liu, J.; Jiang, B. GEO satellite autonomous navigation using X-ray pulsar navigation and GNSS measurements. Int. J. Innov. Comp. Inform. Control 2012, 8, 2965–2977. [Google Scholar]
Figure 1. Model of a satellite in the RSW coordinate system applying the thruster configuration with robotic arms.
Figure 1. Model of a satellite in the RSW coordinate system applying the thruster configuration with robotic arms.
Aerospace 09 00182 g001
Figure 2. Time history of the osculating elements in the one-year simulation for the QFFC.
Figure 2. Time history of the osculating elements in the one-year simulation for the QFFC.
Aerospace 09 00182 g002
Figure 3. Time history of the osculating elements in the one-year simulation for the FOMPC.
Figure 3. Time history of the osculating elements in the one-year simulation for the FOMPC.
Aerospace 09 00182 g003
Figure 4. The comparison of the trajectories in the ( δ λ , δ ϕ ) plane.
Figure 4. The comparison of the trajectories in the ( δ λ , δ ϕ ) plane.
Aerospace 09 00182 g004
Figure 5. The comparison of the trajectories in the ( δ e x , δ e y ) plane.
Figure 5. The comparison of the trajectories in the ( δ e x , δ e y ) plane.
Aerospace 09 00182 g005
Figure 6. The comparison of the trajectories in the ( δ i x , δ i y ) plane.
Figure 6. The comparison of the trajectories in the ( δ i x , δ i y ) plane.
Aerospace 09 00182 g006
Figure 7. The half-year time histories of Δ V NS for the QFFC and the FOMPC.
Figure 7. The half-year time histories of Δ V NS for the QFFC and the FOMPC.
Aerospace 09 00182 g007
Figure 8. The half-year time histories of β ¯ for the QFFC and the FOMPC.
Figure 8. The half-year time histories of β ¯ for the QFFC and the FOMPC.
Aerospace 09 00182 g008
Figure 9. The results of Monte Carlo simulations for β max .
Figure 9. The results of Monte Carlo simulations for β max .
Aerospace 09 00182 g009
Figure 10. The results of Monte Carlo simulations for β min .
Figure 10. The results of Monte Carlo simulations for β min .
Aerospace 09 00182 g010
Figure 11. The results of Monte Carlo simulations for α max .
Figure 11. The results of Monte Carlo simulations for α max .
Aerospace 09 00182 g011
Figure 12. The results of Monte Carlo simulations for α min .
Figure 12. The results of Monte Carlo simulations for α min .
Aerospace 09 00182 g012
Figure 13. The results of Monte Carlo simulations for the prediction horizon N under different levels of constraints.
Figure 13. The results of Monte Carlo simulations for the prediction horizon N under different levels of constraints.
Aerospace 09 00182 g013
Figure 14. The results of Monte Carlo simulations for Δ i max .
Figure 14. The results of Monte Carlo simulations for Δ i max .
Aerospace 09 00182 g014
Figure 15. The results of Monte Carlo simulations for the attitude error of different magnitudes.
Figure 15. The results of Monte Carlo simulations for the attitude error of different magnitudes.
Aerospace 09 00182 g015
Figure 16. The results of Monte Carlo simulations for the determination error of different magnitudes.
Figure 16. The results of Monte Carlo simulations for the determination error of different magnitudes.
Aerospace 09 00182 g016
Figure 17. The trajectories in the ( δ λ , δ ϕ ) plane considering the maximum EDR of the determination.
Figure 17. The trajectories in the ( δ λ , δ ϕ ) plane considering the maximum EDR of the determination.
Aerospace 09 00182 g017
Table 1. The values of β max , β min , α max , α min in different levels of constraints.
Table 1. The values of β max , β min , α max , α min in different levels of constraints.
Level β max β min α max α min
Weak 90 0 90 90
Normal 70 0 60 30
Strong 60 0 45 25
Table 2. The summary of our results and the comparison with those in the previous works.
Table 2. The summary of our results and the comparison with those in the previous works.
V ¯ AN ( m / s / year ) Controlling [ e x , e y ] T
and [ i x , i y ] T
Number of
Thrusters
ON-OFF
Thrusters
Number of
Maneuvers Per
Thruster Per Orbit
QFFC60.2Yes2Yes0.997
FOMPC54.2Yes2Yes0.988
[27]78.6No4Yes-
[30]59.4No6No-
[31]66.8No4Yes2.98
[32]72.3Yes-Yes-
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, C.; Xu, B.; Zhou, W.; Peng, Q. Geostationary Station-Keeping of Electric-Propulsion Satellite Equipped with Robotic Arms. Aerospace 2022, 9, 182. https://doi.org/10.3390/aerospace9040182

AMA Style

Li C, Xu B, Zhou W, Peng Q. Geostationary Station-Keeping of Electric-Propulsion Satellite Equipped with Robotic Arms. Aerospace. 2022; 9(4):182. https://doi.org/10.3390/aerospace9040182

Chicago/Turabian Style

Li, Chengzhang, Bo Xu, Wanmeng Zhou, and Qibo Peng. 2022. "Geostationary Station-Keeping of Electric-Propulsion Satellite Equipped with Robotic Arms" Aerospace 9, no. 4: 182. https://doi.org/10.3390/aerospace9040182

APA Style

Li, C., Xu, B., Zhou, W., & Peng, Q. (2022). Geostationary Station-Keeping of Electric-Propulsion Satellite Equipped with Robotic Arms. Aerospace, 9(4), 182. https://doi.org/10.3390/aerospace9040182

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop