Station-keeping strategy for real translunar libration point orbits using continuous thrust

Aerospace Science and Technology 94 (2019) 105376

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Station-keeping strategy for real translunar libration point orbits using continuous thrust Yi Qi, Anton de Ruiter ∗ Department of Aerospace Engineering, Ryerson University, 350 Victoria Street, Toronto, ON M5B 2K3, Canada

a r t i c l e

i n f o

Article history: Received 27 June 2018 Received in revised form 23 May 2019 Accepted 2 September 2019 Available online 5 September 2019 Keywords: Libration point orbit Station keeping Earth-Moon system Continuous thrust Backstepping technique

a b s t r a c t In this paper, we investigate the station-keeping strategy for translunar libration point orbits (LPOs) using the continuous thrust in the ephemeris model. Ephemeris models with and without the solar radiation pressure (SRP) are proposed for the numerical simulation. Three kinds of translunar LPOs, including halo orbits, vertical Lyapunov orbits, and Lissajous orbits, are used as nominal orbits. A station-keeping strategy based on the backstepping technique is extended to the ephemeris model. Then station keeping under practical constraints, caused by the navigation system and the executive system, is studied in ephemeris models with and without the SRP. Using numerical simulations, influences of the dead-band scheme, navigation errors, the navigation interval time and the area-to-mass ratio on station keeping are discussed, respectively. Furthermore, numerical simulations indicate that if the SRP is taken into account, LPOs obtained in the ephemeris model with the SRP is more preferable as nominal orbits than those without the SRP. More Monte-Carlo simulations testify that our control strategy can be applied to the long-term station keeping for different translunar LPOs. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction Libration point orbits (LPOs) around the Earth-Moon Lagrange Point L 2 , i.e. translunar LPOs, play important roles in both astronautic and scientific fields. According to the Global Exploration Roadmap released by the International Space Exploration Coordination Group, a forum set up by 12 space agencies to advance the global space exploration strategy, a Deep Space Habitat (DSH) near the Moon is required to act as a base station for missions to the Moon or near-Earth asteroids [1]. Translunar LPOs are regarded as the favored choices among international space exploration partners for a DSH. In addition, translunar LPOs could provide continuous communications with both the Earth and the far side of the Moon [1]. Therefore, several space missions of translunar LPOs are proposed or implemented, such as ARTEMIS mission [2] of National Aeronautics and Space Administration (NASA), Chang’e 5-T1 mission of China National Space Administration (CNSA), and CNSA’s Chang’e 4 relay satellite launched in May 2018. However, the dynamics associated with lunar LPOs are unstable, especially the dynamical behavior of translunar LPOs is more complex compared with the other collinear points [3]. Therefore, the station-keeping

*

Corresponding author. E-mail addresses: [email protected] (Y. Qi), [email protected] (A. de Ruiter).

https://doi.org/10.1016/j.ast.2019.105376 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

strategy is more necessary and more challenging for translunar LPOs. There have been numerous researches involving station keeping of LPOs. According to the implementing method of the control strategy, the known station-keeping strategies can be classified into two categories: the impulsive control techniques and the continuous control techniques. The impulsive control techniques are most practical to the present day. According to the Floquet theory, there exist unstable directions in the linear dynamical system of LPOs. An impulsive station-keeping method called Floquet mode approach was proposed to eliminate instability of motion of LPOs by Wiesel and Shelton [4] and Simo et al. [5]. Gomez et al. [6] applied this method to station keeping of translunar LPOs. Hou et al. [7] presented a station-keeping method similar to the Floquet approach to remove the unstable component of the translunar libration point in the real Earth-Moon system. Through this method, the spacecraft can be kept around the translunar libration point. Howell and Pernicka [8] proposed the target point approach (TPA) to minimize a weighted cost function containing fuel costs and state deviations. Lian et al. [9] developed a discrete-time sliding mode control (DSMC) to investigate station keeping for real translunar LPOs in 10 years. Compared with the traditional discrete linear quadratic regulator (DLQR) approach, the parameter tuning of the DSMC is easier. The performance of the DLQR is more sensitive to the values of control parameters than the DSMC strategy. Bai and Junk-

2

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

ins [10] presented the Modified Chebyshev-Picard Iteration (MCPI) method and used it to compute corrective maneuvers at discrete time intervals for station keeping of the halo orbit. Compared with previously published results, their station-keeping approach had a simple control structure that did not require weight tuning and, most importantly, did not need state transition matrix or gradient information computation. Although so far there have been no LPO missions with the continuous control techniques in practice, the use of micro- and nanospacecraft in the study of the Solar System will change this situation in the nearest future [11]. Breakwell et al. [12] applied the linear quadratic regulator (LQR) technique to station keeping for the translunar halo orbit in the circular restricted three-body problem (CRTBP). However, the classical LQR control technique requires a solution to the time-varying Riccati equation [13,14]. To avoid repetitively solving the Riccati equation, a backstepping technique was proposed to reduce the original system to a timeinvariant subsystem for which it is easier to obtain the control law stabilizing this subsystem [15]. Nazari et al. [16] applied the backstepping technique and the dead-band scheme to station keeping of a translunar halo orbit in the Earth-Moon CRTBP. Based on numerical partial derivatives for a reference translunar halo orbit, Ulybyshev [17] proposed an optimization method for the longterm station-keeping of translunar halo orbits. His method is suitable for analysis and simulation of a low-thrust space station in the full ephemeris model. Using a simple linear extended state observer, Narula and Biggs [18] extended a LQR control scheme to improve the station-keeping tracking performance in the presence of disturbances and thruster faults. Based on the nominal multirevolution elliptic halo orbit and the receding horizon control strategy solved by the indirect Radau pseudospectral method, Peng et al. [19] investigated the station-keeping problem in the elliptic restricted three-body problem (ERTBP). Recently, Shirobokov et al. [11] systematically surveyed the known station-keeping techniques for LPOs. According to their survey and discussion, some important recommendations for station keeping were proposed. For example, the nominal orbit of station keeping should be constructed in a high-fidelity model, because a station-keeping method developed in a low-fidelity model of motion normally requires greater costs in high-fidelity models. The uncertainty levels of the orbit determination and maneuver execution processes should be simulated. A series of the Monte Carlo trials is necessary to estimate the efficiency of the station-keeping scheme to give mission designers a clear picture of what to expect from it. In this paper, based on the backstepping technique, we propose a control strategy under practical constraints for this issue. In our control strategy, the navigation interval and the uncertainty from the navigation system and maneuver execution processes are simulated. In addition, the solar radiation pressure (SRP) is also taken into account and analyzed. The influences of practical constraints and nominal orbits will be discussed by Monte-Carlo simulations. Compared with previous references on station keeping of translunar LPOs using continuous thrust [12,16,17], the contributions of this paper are the analysis of the influence of the SRP on station keeping and the extensive numerical investigation for translunar LPOs in the ephemeris model under practical constraints. The structure of this paper can be organized as follows. In Section 2, we introduce the ephemeris model and real translunar LPOs used in this paper. In Section 3, the control method based on the backstepping technique is used and analyzed in the ephemeris model. In Section 4, station keeping of LPOs under practical constraints are discussed by numerical simulations. Finally, the conclusions of this paper are presented in Section 5.

Fig. 1. Moon-centered J2000 inertial coordinates.

2. Background In this section, ephemeris models with and without SRP are introduced. Then, real translunar LPOs are constructed in ephemeris models, which will be used as nominal orbits of the stationkeeping problem in this paper. 2.1. Ephemeris models We use the Moon-centered J2000 inertial coordinates to describe the motion of the spacecraft on translunar LPOs (see Fig. 1). The X axis points to the equinox at the epoch J2000.0, and the X Y plane coincides with the Earth’s equator at the epoch J2000.0. ρ is the position vector of the spacecraft in the Moon-centered J2000 inertial coordinates. P sun and P planet represent position vectors of the Sun and the planet considered in this coordinates, respectively. The position data of celestial bodies in the J2000 inertial coordinates can be obtained from the Jet Propulsion Laboratory (JPL) ephemeris DE430 [20]. If the SRP is not taken into account in the ephemeris model, the equation of motion of the spacecraft in the Moon-centered J2000 inertial coordinates can be expressed by



ρ¨ = −G Mm



 ρ ρ − Pi Pi , − G Mi + G Mi 3 3 ρ  ρ − P i   P i 3 i ∈

(1)

where  = {sun, mercury, venus, earth,..., neptune} is set of the perturbation gravity bodies. In this paper, the Sun and the eight planets are considered. G is the gravitational constant. M m is the mass of the Moon, and M i denotes the corresponding mass of the gravity body in . If the SRP is taken into account, the above equation (1) can be rewritten as



ρ¨ = −G Mm



 ρ ρ − Pi Pi + aS R P , − G Mi + G Mi ρ 3 ρ − P i 3  P i 3 i ∈ (2)

where a S R P is the SRP acceleration and can be obtained by [21]

a S R P = ηβ

G M sun

ρ − P sun 2

(cos α )2 n,

(3)

where n is the unit normal vector of the solar array or solar sail (see Fig. 1).

n=

ρ − P sun cos α + t sin α , ρ − P sun 

(4)

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

where t is the unit vector perpendicular to the vector from the Sun to the spacecraft (see Fig. 1). α is the angle between the normal vector n and the vector from the Sun to the spacecraft, and α ∈ [0o , 90o ]. η in Eq. (3) is the reflectivity of the solar array or solar sail, typically around 0.9 for a real spacecraft, with an ideal maximum value of 1. β is the solar array or solar sail lightness parameter defined as the ratio of the SRP force to the solar gravitational attraction [21]

β=

Lsλ 2π cG M sun

≈ 1.5441 × 10−3

A m

(5)

,

where L s = 3.86 × 1026 W is the solar luminosity, c = 2.99792 × 108 m/s is the speed of light, and λ = A /m is the ratio of surface area to mass. In the numerical calculation of this paper, we use parameters of the SOHO spacecraft to calculate λ, where m = 1850 kg, and A = 9.5 m × 3.65 m = 34.675 m2 [22], so λ = 0.0187 m2 /kg. In practical missions, α should be as small as possible to take full advantage of the solar power or the SRP force. If we assume that α = 0o , based on Eqs. (2)∼(5), we can get



ρ¨ = − G Mm

 ρ ρ − Pi Pi − G Mi + G Mi 3 ρ  ρ − P i 3  P i 3 i ∈

+ ηβ G M sun = − G Mm

ρ − P sun ρ − P sun 3

ρ ρ 3

−

 i ∈

G Mi



(6)



ρ − Pi Pi , + G Mi ρ − P i 3  P i 3

where M i = (1 − ηβ) M i if i = sun; otherwise, M i = M i . We find that the right hand of Eq. (6) is quite similar to that of Eq. (1), so they are denoted by f s (t , ρ ) and f (t , ρ ), respectively. 2.2. Real translunar libration point orbits In this paper, three kinds of translunar LPOs, including halo orbit, Lissajous orbit and vertical Lyapunov orbit, are considered. Ephemeris models with and without the SRP acceleration are both taken into account. The LPOs obtained in this subsection will be used as nominal orbits of station keeping in following sections. Numerical construction methods of LPOs have been well investigated by many researchers [23–25], and they are not our emphasis of this paper. Hence, in this subsection, the construction procedure is briefly summarized. Readers can refer to corresponding literature for details. The symmetric periodic LPOs, such as halo orbits and vertical Lyapunov orbits, can be constructed in the ephemeris model using four steps as follows [24],

• Step 1: Generate the third order approximation of LPOs using the expressions presented by Richardson [26];

• Step 2: Use the third order approximation obtained in Step 1 as the initial guess, and construct symmetric periodic LPOs in the Earth-Moon CRTBP using the single-shooting differential correction proposed by Howell [23]; • Step 3: Use the LPOs obtained in Step 2 as the initial guess, and construct corresponding LPOs in the ephemeris model by the multiple-shooting differential correction [27,24]; • Step 4: Prune the first and last orbits of the LPOs obtained in Step 3 in order to observe the structure of the resulting quasi-periodic orbit [24]. For the quasi-periodic orbit, such as the Lissajous orbit, it can also be obtained by the above process except Step 2, so the initial

3

guess of the Lissajous orbit in Step 3 directly adopt the result of Step 1. In numerical computations of Step 3, Eq. (1) is used to calculate the ephemeris model without the SRP acceleration, and for ephemeris model with the SRP acceleration, we assume η = 0.9 and use Eq. (6), where the angle α = 0o . The numerical integrator in the multiple shooting approach is MATLAB’s function ode113, which is a variable-step, variable-order Adams-Bashforth-Moulton PECE solver of orders 1 to 13. The absolute and relative error tolerances are both set as 10−9 . Fig. 2 shows halo orbits, vertical Lyapunov orbits, and Lissajous orbits obtained by the above process, where red and blue lines denote the results of ephemeris models with and without the SRP acceleration, respectively. Those orbits are described in the dimensionless Earth-Moon rotating frame where the origin is located at the Earth-Moon barycenter, the xem axis points to the Moon, the xem y em plane coincides with the Moon’s orbit, and the length unit is the instantaneous distance between the Earth and the Moon. Hence, the position of the Moon is (1 − μ, 0, 0) in the dimensionless Earth-Moon rotating frame, where μ = M m /( M earth + M m ) = 0.01215058 is the mass ratio of Earth-Moon system. The libration point L 2 is located at (1.15568248, 0, 0). The data of LPOs in Fig. 2 are listed in Table 1, where A z is the maximum amplitude along zem . The durations of LPOs in ephemeris models are about 2 years (1 year = 365 days). These LPOs in Fig. 2 are used as nominal orbits of the stationkeeping problem in this paper. As we can see from Fig. 2, shapes and positions of LPOs with the SRP are quite similar to the corresponding orbits without the SRP. To compare the influence of the SRP, Fig. 3 indicates the position and velocity differences of LPOs with and without the SRP at different time, where black, blue, and red lines denote results of halo orbits, vertical Lyapunov orbits, and Lissajous orbits, respectively. Numerical calculation indicates that averages of the position differences of the halo orbits, vertical Lyapunov orbits, and Lissajous orbits are 15.4898, 13.0902, and 12.5869 km, respectively; averages of the velocity differences of the halo orbits, vertical Lyapunov orbits, and Lissajous orbits are 0.0572, 0.0536, and 0.0510 m/s, respectively. In the station-keeping problem, these differences cannot be ignored, i.e., the influence of the SRP should be taken into account in practical LPO missions. As seen in Fig. 3, the SRP has an obvious influence on LPOs. According to Eqs. (3) and (5), the SRP acceleration is determined by the area-to-mass ratio of the spacecraft λ. To show the influence of λ on LPOs, we choose different values of λ to construct halo orbits using the above approach. Fig. 4 shows the position and velocity differences of nominal halo orbits with different values of λ compared to the halo orbit without the SRP, where red, black, blue, and green lines denote results of λ = 0.01, 0.0187, 0.03, and 0.04 m2 /kg, respectively. From the figure, we can observe that with the increase of λ, both position and velocity differences increase. Numerical computations indicate that averages of the position differences of the halo orbits with λ = 0.01, 0.0187, 0.03, and 0.04 m2 /kg, are 8.2640, 15.4898, 24.7934, and 33.0588 km, respectively; averages of the velocity differences of the halo orbits are 0.0305, 0.0572, 0.0916 and 0.1221 m/s, respectively. 3. Backstepping technique In this section, the control strategy for the continuous thrust based on the backstepping technique is extended to the stationkeeping problem of LPOs in the ephemeris model. 3.1. Control strategy The controlled equation of motion of the spacecraft in the ephemeris model can be obtained by

4

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

Fig. 2. Three kinds of LPOs: (a) halo orbits, (b) vertical Lyapunov orbits, (c) Lissajous orbits, in the dimensionless Earth-Moon rotating frame, where blue lines denote LPOs in the real ephemeris model without the SRP, and red lines denote LPOs in the real ephemeris model with the SRP. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Table 1 Data of LPOs in Fig. 2. LPOs

Start time (UTC)

End time (UTC)

A z (km)

Halo without SRP Halo with SRP Vertical Lyapunov without SRP Vertical Lyapunov with SRP Lissajous without SRP Lissajous with SRP

2018 2018 2018 2018 2018 2018

2020 2020 2020 2020 2020 2020

4.3362 × 104 4.3362 × 104 8.4259 × 104 8.4257 × 104 2.7079 × 104 2.7078 × 104

JUN JUN JUN JUN JUN JUN

14 14 12 12 13 13

14:25:45 14:25:45 23:39:43 23:39:43 08:43:23 08:43:23

JUN JUN JUN JUN JUN JUN

13 13 11 11 12 12

14:28:24 14:28:24 23:40:55 23:40:55 08:43:44 08:43:44

Fig. 3. Differences of nominal orbits with and without the SRP, where black, blue, and red lines denote results of halo orbits, vertical Lyapunov orbits, and Lissajous orbits, respectively.

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

5

Fig. 4. Differences of halo orbits with different area-to-mass ratio λ compared to the halo orbit without the SRP.

x˙ = f (s) (t , x) + Bu (t ),

(7)

where x is the phase state of the spacecraft in the Moon-centered J2000 inertial coordinates. B = [ O 3×3 I 3×3 ] T , where O 3×3 and I 3×3 are 3 × 3 zero matrix and identity matrix, respectively. u ∈ R3 is the control input or the thrust acceleration. The state deviation with respect to the nominal orbit is x = x − x N , where x N is the phase state of the nominal orbit, so we can get x˙ N = f (s) (t , x N ). In this way, we can derive the linear dynamical equation of the state deviation as follows

x˙ = A (t , x N )x + Bu (t )    O 3×3 I 3×3 ρ + Bu (t ) = ∂ f (s) /∂ ρ O 3×3 v    O 3×3 I 3×3 ρ + Bu (t ), = A 21 (t , ρ ) O 3×3 v

x˙ 1 = x2 , x˙ 2 = A 21 (t )x1 + u .

(8)

(10)

Hence, taking the derivative of Eq. (9) with respect to t, we can get

z˙ 1 = S z 1 + z 2 , z˙ 2 = ( A 21 (t ) − S 2 ) z 1 − S z 2 + u .

(11)

Let u = −( A 21 − S 2 ) z 1 + u 1 , then, the time-varying term in the second equation of Eq. (11) can be removed

z˙ 2 = − S z 2 + u 1 (t ).

where A is the Jacobi matrix. For the ephemeris model without the SRP, A 21 (t , ρ ) = ∂ f /∂ ρ ; for the ephemeris model with the SRP, A 21 (t , ρ ) = ∂ f s /∂ ρ . Since the phase state of the nominal orbit x N (t ) is given, A (t , x N ) = A (t ) only depends on time. The classical optimal linear quadratic regulator (LQR) control technique to minimize the quadratic cost function

(12)

For the time-invariant subsystem (12), we assume that u 1 (t ) = −kz 2 (t ), where k should satisfy the minimization of the quadratic cost function

∞ J = [ z 2T (t ) Q z 2 (t ) + u 1T (t ) Ru 1 (t )]dt ,

(13)

0

t f [xT (t ) Q x(t ) + u T (t ) Ru (t )]dt ,

J=

According to Eq. (8), we can obtain

t0

requires a solution to the time-varying Riccati equation. To avoid repetitively solving the Riccati equation, we use the backstepping technique proposed by Nazari et al. [16] to reduce the original system into a time-invariant subsystem for which it is easier to obtain the control law stabilizing this subsystem. Based on the result of the subsystem, we can derive the control law that stabilizes the original system by working backward. According to references [15, 16], the backstepping transformation is expressed by

z 1 = x1 , z 2 = x2 − S x1 ,

(9)

where x1 = ρ − ρ N and x2 = v − v N are position and velocity deviations with respect to the nominal orbit, respectively. S is a constant negative definite matrix.

where Q and R are constant weighting matrices. Q ∈ R3×3 is a positive semi-definite matrix, and R ∈ R3×3 is a positive definite matrix. To calculate k, we need to solve an algebraic Riccati equation (ARE) for the infinite horizon problem

S T P + P S − Q + P R −1 P = O 3×3 .

(14)

Then, the control gain of the subsystem (12), k, can be obtained by k = R −1 P . Therefore, in the numerical computation of k, we just need to solve an ARE once. For the original system (7), the control input can be expressed as

u = − ( A 21 − S 2 ) z 1 − kz 2 (t )

= − ( A 21 − S 2 + k S )x1 − kx2 = − K x,

(15)

where A 21 can be obtained from the nominal orbit (see Eq. (8)).

6

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

3.2. Idealistic simulation In this subsection, we use the idealistic simulation to choose values of parameters in the backstepping technique, where practical errors, such as navigation errors or the maneuver error, are not taken into consideration, except initial insertion errors. By ignoring all the other disturbing factors, one has a neat picture of the performance of the controller [16]. Before the simulation, we analyze how to choose the appropriate values of parameters in the backstepping technique, including S , Q , and R. Firstly, we estimate the magnitude of S . In reference [3], the initial insertion error includes the position error of 1 km and the velocity error of 10−5 km/s (= 1 cm/s). Hence, we assume that the magnitude of x2 is about 10−5 km/s, and the magnitude of x1 is about 1 km. Since z 1 and z 2 in Eq. (9) are the equivalent parameters of x1 and x2 , respectively, based on the second equation of Eq. (9), we estimate that the magnitude of S should be smaller than 10−5 . Then, we discuss the relationship of Q and R. Based on Eq. (13), the terms of z 2T Q z 2 and u 1T R u 1 should have similar magnitudes. As previous discussion, the magnitude of z 2 is same as x2 , about 10−5 km/s (= 1 cm/s). For the usual low-thrust engine with thrust of 1 mN and a common spacecraft with mass of 2000 kg, the magnitude of the thrust acceleration u is about 10−10 km/s2 . Then, we estimate that (10−5 )2  Q  ≈ (10−10 )2  R . According to the above discussion, we set S = −10−6 I 3×3 , Q = 10−4 I 3×3 , and R = 106 I 3×3 . In numerical simulations of this paper, initial insertion errors are assumed to be white noise processes and have standard normal distributions with zero means and standard deviations of σ X ,Y , Z = 100 km and σ X˙ ,Y˙ , Z˙ = 1 cm/s. The result of an idealistic simulation is displayed in Fig. 5, where we adopt the ephemeris model without the SRP. The simulation time is one year, and the nominal orbit is the halo orbit without the SRP in Fig. 2. In this paper, the numerical integrator in the simulation of station keeping is MATLAB’s function ode45, which is based on an explicit Runge-Kutta (4,5) formula, the DormandPrince pair. The absolute and relative error tolerances are both set as 10−6 . In this paper, d(t ) = [x1 (t ) x2 (t ) x3 (t )] T denotes the position error vector of the spacecraft to the nominal orbit, and the position deviation d(t ) = d(t ). As we can see from Fig. 5(a), position errors decrease to about zero within about 50 days. Similar to previous reference [16], we also adopt the 2% settling time envelope of the position error to verify the rapidity of the control strategy. In [16], the 2% settling time is about 0.27T , where T is the simulation time. Numerical calculation indicates that in Fig. 5(a), the time that it takes for the displacement between the spacecraft and the nominal orbit to approach 2% of the initial state deviation is equal to 46.3550 days (<0.2 years). The result of the control input vector u (t ) = [u 1 (t ) u 2 (t ) u 3 (t )] T is illustrated in Fig. 5(b), where the magnitude of the thrust acceleration u (t ) = u (t ). Apparently, u also decreases to about zero within about50 days. The t total fuel cost of station keeping is defined as  v = 0 u (t )dt, and the  v in Fig. 5 is 5.3138 m/s. 3.3. Dead-band scheme To avoid thrusters working continuously for long-term maneuvers, the authors in [16] proposed a dead-band scheme for the station-keeping issue using the continuous thrust. The dead-band control scheme is quite simple, and is illustrated in Fig. 6. It assumes that thrusters shut down when the distance of the spacecraft from the nominal orbit decreases to a lower boundary from a farther position (see t 1 in Fig. 6), and then thrusters start up when

Fig. 5. An idealistic simulation of station keeping for the halo orbit without the SRP: (a) position errors versus time t, (b) control input versus t.

Fig. 6. Dead-band control schematics.

the distance from the nominal orbit increases to a upper boundary from a closer position (see t 2 in Fig. 6). Obviously, the interval time between t 1 and t 2 is the idle time of thrusters. If the upper and lower boundaries are set as 20 km and 5 km, respectively, Fig. 7 shows the result of a numerical example of the backstepping technique with the dead-band scheme, where the green bars denote idle segments of dead-band scheme. The ephemeris model without the SRP is considered. To compare the influence of the dead-band scheme, the nominal orbit of the simulation is also the halo orbit without the SRP in Fig. 2, and the simulation time is still one year, similar to the simulation in the previous subsection. By looking at the figure, we find that the po-

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

7

Fig. 7. Control strategy with dead-band scheme: (a) position deviation to the nominal orbit, (b) control input.

sition deviation d decreases to about 5 km within 50 days, then fluctuates between 5 km and about 25 km when the dead-band scheme works. Numerical computation indicates that the fuel cost  v = 15.3677 m/s in Fig. 7(b) is larger than that in Fig. 5(b), which accords with the observation of two figures. At the expense of a larger fuel cost, the total idle time in this figure is T idle = 81.7996 days. Next, we analyze influences of upper and lower boundaries of the dead-band scheme on the fuel cost  v and idle time T idle using Monte-Carlo simulations. In the Monte-Carlo simulations, for each case with the given upper and lower boundaries, 20 sample points with different initial insertion errors are chosen. As mentioned before, initial insertion errors are white noise processes and have standard normal distributions with zero means and standard deviations of σ X ,Y , Z = 100 km and σ X˙ ,Y˙ , Z˙ = 1 cm/s. Similar to the idealistic simulation in the previous subsection, only the initial insertion error is taken into account. The nominal orbit of the simulation is also the halo orbit without the SRP in Fig. 2, and the simulation time is still one year. The ephemeris model without the SRP is taken into account. We assume that the lower boundary is fixed and set as 5 km, but the upper boundary is varying. Performances of different upper boundaries are displayed in Fig. 8, where black points denote average values of corresponding Monte-Carlo simulations. From these figures, we find that the average fuel cost  v increases with an increase in upper boundaries, but the average idle time T idle does not have a such well-behaved rule. In this paper, we postulate that a good station-keeping strategy should have a smaller fuel cost  v and a larger idle time T idle . Based on this requirement, we postulate that the upper boundary should be 10∼20 km.

Fig. 8. Performances of different upper boundaries: (a) fuel cost  v, (b) idle time T idle .

If we assume that the upper boundary is fixed and set as 20 km, but the lower boundary is varying. Performances of different lower boundaries are displayed in Fig. 9, where black points denote average values of corresponding Monte-Carlo simulations. It should be noted that the situation of 0 km actually corresponds to the control strategy without the dead-band scheme. As we can see, the average fuel cost  v gradually increases with an increase in lower boundaries. The control strategy without a dead-band scheme apparently has no idle times, but for the dead-band strategy, T idle s are distributed in the range of 60∼80 days, where the maximum T idle belongs to lower boundary of 5 km. To show the influence of the upper and lower boundaries in detail, we chose a fixed initial insertion error: [x1 (t ), x2 (t ), x3 (t )] = [32.64, 68.83, −44.49] km and [x˙ 1 (t ), x˙ 2 (t ), x˙ 3 (t )] = [−0.809, 0.179, −0.466] cm/s. Fig. 10 displays the distributions of the fuel cost  v and idle time T idle with different upper and lower boundaries under this initial insertion error. As we can see from Fig. 10, the fuel cost  v generally increases with the increase of upper and lower boundaries, but the distribution of the idle time T idle is more complex, which is similar to the results in Figs. 8 and 9. We observe that when the lower boundary is smaller than 1 km, the fuel cost is small (about 5 m/s), but the corresponding idle time is lower than 30 days. Based on the results in Figs. 8∼10, we find that the dead-band scheme with an upper boundary of 20 km and a lower boundary of 5 km is a good option when the fuel cost and the idle time are both taken into account. Therefore, we use these two values in the dead-band scheme of the following simulations, which will not be repeated again.

8

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

4.1. Control strategy with practical constraints In this subsection, the control strategy with practical constraints is proposed in the ephemeris model. In actual LPO missions, the navigation system does not provide updates in continuous time, so the navigation interval should be taken into account. We use t i and  T n to denote the navigation time and the navigation interval time, respectively. Once the state error with respect to the nominal orbit x(t i ) = xni is obtained from the navigation system, then based on Eq. (8), in the navigation interval [t i , t i +  T n ], the state error can be obtained by forward numerical integration of the linear dynamical system as follows,

x˙ = A (t , xN )x + Bu (t ), x(t i ) = xni .

(16)

Note that x is obtained from the linear dynamical system (16) rather than the real nonlinear ephemeris model, so it actually is the pseudo state error rather than the real state error x − x N . As mentioned before, A (t , x N ) = A (t ) only depends on time, then in the numerical integration of above equation (16), we do not need any other information from the navigation system except xni . In addition, since A (t ) could be obtained a prior, the computational cost of the numerical integration of Eq. (16) is lower than that of the real ephemeris model, i.e. the former is more computationally implementable on a flight processor. Based on Eq. (15), the idealistic control input in the navigation interval [t i , t i +  T n ] can be expressed as

u¯ (t ) = − K (t , x N )x, Fig. 9. Performances of different lower boundaries: (a) fuel cost  v, (b) idle time T idle .

4. Station keeping with practical constraints In this section, station keeping of LPOs under practical constraints is investigated. Practical constraints considered here come from the navigation system and the executive system. Constraints of the navigation system include navigation errors and the navigation interval. Constraints of the executive system include the maneuver error and the engine limitation. Ephemeris models with and without the SRP are both used and studied in numerical simulations.

(17)

where x is the pseudo state error obtained from Eq. (16). It should be noted that both the engine limitation and the maneuver error are not considered in Eq. (17). We assume that the engine limitation requires that the magnitude of the thrust acceleration performed should be constrained in the range of [umin , umax ]. If u¯ = u¯  < umin , u¯ = 0. Hence, u¯ (t ) can be rewritten as ∗

u (t ) =

⎧ ⎨ ⎩

0,

if

u¯ < umin ;

umax u¯ , u¯

if

u¯ > umax ;

u¯ ,

if

(18)

u¯ ∈ [umin , umax ].

According to [9], the maneuver error is due to the imperfection of the manufacturing or installation of the thrusters, which leads to

Fig. 10. Performances of different upper and lower boundaries: (a) fuel cost  v, (b) idle time T idle .

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

Table 2 Parameters for practical constraints. Parameters

Table 3 Comparison of station keeping under practical constraints with and without deadband strategy.

Values

Initial insertion errors

σ X ,Y , Z σ X˙ ,Y˙ , Z˙

100 km 1 cm/s

Navigation errors

σ Xn ,Y n , Z n σ X˙ n ,Y˙ n , Z˙ n

1 km 1 cm/s

Maneuver errors

σu X , u Y , u Z

2%

Navigation interval T n

2 days

Engine limitation umax umin

5 × 10−4 m/s2 1 × 10−7 m/s2

inaccuracy in both the magnitude and the direction of u (t ). We assume that the maneuver error is denoted by εu i , i = X , Y , Z . Then, the control input executed in Eq. (16) is



u (t ) =

0,

if

E m u∗ ,

if

u ∗ < umin ; u ∗  umin ,

(19)

where E m = diag [1 + εu X , 1 + εu Y , 1 + εu Z ]. If the dead-band scheme is considered in station keeping, u = 0 when the deadband condition occurs. The navigation error is denoted by εn = [ε X εY ε Z ε X˙ εY˙ ε Z˙ ] T , then, the actual state error obtained from the navigation system at the time t i can be expressed as,

xni = xreal (t i ) + εn ,

9

(20)

where xreal (t i ) is the real state error without errors. Using Eqs. (16)-(20), in the navigation interval [t i , t i +  T n ], the control strategy with practical constraints is established. In the following numerical simulations, we assume that both the navigation error and maneuver error are white noise processes and have standard normal distributions with zero means. Their standard deviations and other values of parameters are listed in Table 2. In addition, as mentioned before, initial insertion errors are white noise processes and have standard normal distributions with zero means. The parameters of the navigation system come from the navigation uncertainties of ARTEMIS mission [28]. 4.2. Numerical simulations without SRP In this subsection, the real ephemeris model without the SRP, i.e. Eq. (1), is used in the station-keeping problem, and LPOs without the SRP in Fig. 2 are used as nominal orbits. Firstly, we investigate the influence of the dead-band scheme on station keeping with practical constraints by the Monte-Carlo simulations. In the Monte-Carlo simulations, for each simulation with or without dead-band strategy, 100 sample points with different initial insertion errors and practical constraint parameters are chosen. As mentioned before, the initial insertion error, navigation error and maneuver error are white noise processes and have standard normal distributions with zero means. Their standard deviations are listed in Table 2. The nominal orbit of the simulation is the halo orbit, and the simulation time is one year. The fuel cost  v and the total idle time T idle are two important indexes considbe noted that in the following research,  v ered here. It should t t is defined as 0 u ∗ (t )dt rather than 0 u (t )dt, because according to the definition of u in Eq. (19), u (t ) randomly fluctuates around u ∗ (t ) within a narrow range. For convenience, we use u ∗ (t ) instead of u (t ) to compute  v and show the control input. For the

Monte-Carlo results

Dead-band

Non dead-band

Avg. v (m/s) Max. v (m/s) Min. v (m/s) Std. v (m/s) Avg.T idle (days) Max.T idle (days) Min.T idle (days) Std.T idle (days)

18.7485 25.7635 13.8800 3.2189 50.7733 62.8617 35.6931 8.02412

14.8686 19.8153 10.5537 2.5755 22.1662 32.1885 10.5273 5.5177

control strategy with a dead-band scheme, the total idle time T idle includes two parts: the idle time of the dead-band scheme and the idle time of u¯ < umin , but for the control strategy without a dead-band scheme, there is only the latter part. The Monte-Carlo results are displayed in Table 3, where Avg., Max., Min., and Std. are abbreviations of the words average, maximum, minimum, and standard deviation, respectively. As we can see from average values in the table, using just about 4 m/s more fuel cost, the control strategy with a dead-band scheme can have more than twice idle time than that without a dead-band scheme. Therefore, synthesizing two aspects of the fuel cost and the idle time, we postulate that dead-band schemes are better than non dead-band schemes for the station-keeping strategy. In the following research of station keeping under practical constraints, the dead-band scheme is the default option for the control strategy. Fig. 11 shows the detail of station keeping with practical constraints and a dead-band scheme. The nominal orbit of the simulation is the halo orbit, and the simulation time is one year. Numerical calculation indicates that fuel cost  v = 15.1010 m/s, and the total idle time T idle = 50.0766 days. The green and pink bars in Fig. 11(a)∼(c) are idle segments causing by the dead-band scheme and umin , respectively. dreal denotes the real position deviation with respect to the nominal orbit, and d pse is the pseudo position deviation obtained by Eq. (16). As we can see from Fig. 11(a), dreal firstly decreases to about 5 km within 50 days, then fluctuates in the range of 3∼35 km. It usually increases rapidly in dead-band idle segments, and decreases in other segments. The curve of d pse fluctuates around that of dreal in a small range. From Fig. 11(b), the detail of the change of the position deviation is displayed. At the navigation time (see black points in that figure), the state of spacecraft with navigation errors is obtained by the navigation system. Then, during the navigation interval, the control strategy with practical constraints proposed in the previous subsection is conducted. As we can see, d pse obtained from the linear dynamical system gradually leaves from dreal with time. But at the next navigation time, the deflected state of the spacecraft will be partially corrected by the navigation system. In addition, the dead-band scheme works based on the pseudo result d pse . From Fig. 11(a) and (b), we find that idle times of the dead-band scheme (about 7 days) are larger than those of u¯ < umin , although the latter are more frequent than the former. From Fig. 11(c), we observe that control input u ∗ decreases from about 8 × 10−6 m/s2 to umin within 50 days, then it fluctuates in the range of about 0∼3 × 10−6 m/s2 . The actual controlled halo orbit is displayed in Fig. 11(d), where red and green lines denote idle and thrust segments, respectively. Another simulation example is illustrated in Fig. 12, where the nominal orbit of the simulation is the Lissajous orbit, and the simulation time is also one year. Its fuel cost  v = 13.9270 m/s, and its total idle time T idle = 58.3300 days. By looking at the figure, station keeping of the Lissajous orbit has similar conclusions with that of the halo orbit in Fig. 11. Furthermore, we can observe idle segments of the controlled Lissajous orbit more clearly in Fig. 12(c).

10

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

Fig. 11. Control strategy with dead-band scheme for the halo orbit: (a) position deviation to the nominal orbit, (b) detail of the position deviation, (c) control input, (d) controlled halo orbit in the Earth-Moon rotating frame, where red and green lines denote idle and thrust segments, respectively.

Fig. 12. Control strategy with dead-band scheme for the Lissajous orbit: (a) position deviation to the nominal orbit, (b) control input, (c) controlled Lissajous orbit in the Earth-Moon rotating frame, where red and green lines denote idle and thrust segments, respectively.

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

11

Fig. 13. Control strategy without navigation errors for the halo orbit: (a) position deviation to the nominal orbit, (b) control input, (c) controlled halo orbit in the Earth-Moon rotating frame, where red and green lines denote idle and thrust segments, respectively.

Next, we analyze the influence of navigation errors on station keeping. Fig. 13 shows the results of station keeping of the halo orbit without navigation errors. As we can see from Fig. 13(a), the curve of d pse is very close to that of dreal . Therefore, we conclude that although the control strategy proposed in the previous subsection is based on the linear dynamical system (16) and the navigation result, the pseudo state of the spacecraft can well approximate the real state of the nonlinear dynamical system if the navigation result is accurate. In other words, the result of station keeping with other practical constraints except navigation errors actually is similar to that without practical constraints in Sect. 3.3, which can be verified by comparing Figs. 7 and 13. Then, we can use the Monte-Carlo results in Sect. 3.3 as results of station keeping for the halo orbit without navigation errors. In Fig. 8, when the upper and lower boundaries are set as 20 km and 5 km, we calculate the average  v = 18.3207 m/s, and the average T idle = 78.8021 days. Based on Table 3, for the control strategy with practical constraints and the dead-band scheme, the average  v = 18.7485 m/s, and the average T idle = 50.7733 days. Hence, we postulate that navigation errors in station keeping with practical constraints can result in less idle time T idle , but have no obvious influence on the fuel cost  v. For example, we can observe that idle times of the deadband scheme in Fig. 13 are larger than those in Fig. 11. The idle interval time is about 12 days in Fig. 13. The actual controlled halo orbit is displayed in Fig. 13(c). Numerical computation indicates that the fuel cost  v = 19.3435 m/s, and the total idle time of thruster T idle = 73.9756 days. 4.3. Numerical simulations with SRP In this subsection, the real ephemeris model with the SRP acceleration, Eq. (2), is used in the station-keeping problem under practical constraints. As mentioned before, in practical missions,

the value of α is small to take full advantage of the solar power or the SRP force. In numerical simulations of this subsection, we assume that η = 0.9, α has a uniform distribution in the range of [0o , 30o ], and the unit vector t in Eq. (4) is uniformly distributed at random in the plane perpendicular to the vector from the Sun to the spacecraft. Both α and t change with time in a statistically independent way. Firstly, we compare results of station keeping using nominal orbits with and without the SRP. For example, we choose halo orbits with and without the SRP in Fig. 2 as nominal orbits, respectively. The simulation time is one year. Fig. 14 shows their control strategies: (a) and (b) are the position deviation and the control input using the nominal halo orbit without the SRP; (c) and (d) are results using the nominal halo orbit with the SRP. As we can see from Fig. 14(a), the curve of dreal decreases to about 25 km within 40 days, and then it is maintained in a range around 25 km near one year. d pse still fluctuates around dreal in a small range. Since the value of d pse is larger than the lower boundary of the dead-band scheme, 5 km, the dead-band condition does not occur during one year. Its narrow idle segments all come from the condition of u¯ < umin . However, according to Fig. 14(c), dreal and d pse using the nominal orbit with the SRP are distributed in a lower range, 5∼25 km. Dead-band idle segments exist, so its total idle time T idle is larger than that using the nominal orbit without the SRP. In addition, comparing Fig. 14(b) and (d), we observe that control strategy corresponding to the halo orbit without the SRP need a larger control input u ∗ and more frequent thrust segments to track the nominal orbit. Numerical computation verifies our conclusions: the fuel cost  v and the total idle time of thrusters T idle for Fig. 14(a) and (b) are 31.2685 m/s and 0.6841 days, respectively; the fuel cost  v and the total idle time T idle for Fig. 14(c) and (d) are 19.6650 m/s and 34.5071 days, respectively.

12

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

Fig. 14. Control strategies with the SRP: (a) and (b) are the position deviation and control input using the nominal halo orbit without SRP; (c) and (d) are the position deviation and control input using the nominal halo orbit with the SRP. Table 4 Comparison of station keeping using different nominal halo orbits.

Table 5 Comparison of station keeping using different nominal LPOs with the SRP.

Monte-Carlo results

Without SRP

With SRP

Monte-Carlo results

Halo

Vertical Lyapunov

Lissajous

Avg. v (m/s) Max. v (m/s) Min. v (m/s) Std. v (m/s) Avg.T idle (days) Max.T idle (days) Min.T idle (days) Std.T idle (days)

33.0038 39.4129 26.9921 3.0622 0.8583 1.7442 0.2423 0.4976

18.6002 23.6432 14.8285 2.9872 45.1722 69.9073 24.8129 12.9095

Avg. v (m/s) Max. v (m/s) Min. v (m/s) Std. v (m/s) Avg.T idle (days) Max.T idle (days) Min.T idle (days) Std.T idle (days)

18.6002 23.6432 14.8285 2.9872 45.1722 69.9073 24.8129 12.9095

16.4707 19.6238 13.4730 1.6721 43.0792 58.5108 31.3389 8.5963

17.8681 22.9624 14.1457 2.7003 43.4196 59.6923 26.2617 8.8751

More numerical simulation results using different nominal halo orbits are listed in Table 4. In the Monte-Carlo simulations, for the given nominal halo orbits, 100 sample points with different α , initial insertion errors and practical constraint parameters (see Table 2) are chosen. From Table 4, we find that performances, including  v and T idle , of station keeping of the nominal halo orbit with the SRP are obviously better than those without the SRP. Therefore, we conclude that if the SRP is taken into consideration in the station-keeping problem, nominal LPOs obtained in a more precise ephemeris model with the SRP are more preferable than those without the SRP. This conclusion verifies the importance of modeling accurate nominal orbits, which accords with the result in [29]. In addition, it should be noted that nominal LPOs with the SRP in Fig. 2 are obtained in Eq. (6), where α is fixed and set as 0o , but in numerical simulations of this subsection, Eq. (2) is used, where α is uniformly distributed at random in the range of [0o , 30o ]. Although ephemeris models used in nominal orbits and station keeping are not identical, based on the result in Table 4, we postulate that those nominal LPOs with the SRP in Fig. 2 are acceptable and valid. Comparing the Monte-Carlo results in the first column of Table 3 and the second column of Table 4, we observe that the addition of the SRP in station keeping has no obvious influence on the fuel cost, but reduces the total idle time from

about 50 days to about 45 days. We postulate that this reduction is attributed nonidentical SRP ephemeris models used in nominal orbits and station keeping. Results of the Monte-Carlo simulations using nominal LPOs with the SRP in Fig. 2 are compared in Table 5. For each simulation with the given nominal LPOs with the SRP, 100 sample points with different α , initial insertion errors and practical constraint parameters (see Table 2) are chosen. The simulation time is one year. From the data in this table, we find that whether in the fuel cost  v or in the total idle time T idle , there is no obvious difference among different LPOs. Therefore, the control strategy under practical constraints we proposed can be generally applied to translunar LPOs in the ephemeris model with the SRP. In Sect. 4.1, the navigation interval time  T n is set as 2 days, i.e., the state of the spacecraft is obtained using nominal tracking arcs of one contact every other day, similar to ARTEMIS Mission. Next, we use the Monte-Carlo simulation to analyze the influence of  T n on the performance of station keeping under practical constraints. In the numerical simulation, the halo orbit with the SRP in Fig. 2 is chosen as the nominal orbit, and the simulation time is one year. Fig. 15 shows the performance of different navigation interval time  T n . From the figure, we find that the fuel cost  v gradually increases with an increase in  T n , but the total idle time  T n decreases with

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

13

Table 6 Comparison of station keeping of halo orbits with different area-to-mass ratio λ. Monte-Carlo results

λ (m2 /kg) 0.01

0.0187

0.03

0.04

Avg. v (m/s) Max. v (m/s) Min. v (m/s) Std. v (m/s) Avg.T idle (days) Max.T idle (days) Min.T idle (days) Std.T idle (days)

19.3076 24.9872 13.8675 3.0898 48.3907 60.9854 33.5295 8.1621

18.6002 23.6432 14.8285 2.9872 45.1722 69.9073 24.8129 12.9095

20.0679 25.2635 14.7581 3.2569 35.0331 46.6921 19.0494 7.2063

19.5772 25.9388 15.2432 2.8494 28.1886 57.1853 9.7959 12.9395

Table 7 Results of station keeping using different nominal LPOs with the SRP for 2 years.

Fig. 15. Performances of different navigation interval time  T n : (a) fuel cost  v, (b) idle time T idle .

the increase of  T n . Apparently, we can conclude that for  T n in station keeping, the shorter the better. However, in practical missions,  T n usually is limited by the ability of the trajectory determination system in ground stations, and the shorter  T n is more challenging for the hardware facility. As mentioned before, the area-to-mass ratio λ can significantly influence halo orbits in the ephemeris model (see Fig. 4), so we investigate the results of station keeping for halo orbits with different λ. We use the halo orbits with different λ as nominal halo orbits, then the Monte-Carlo simulations are implemented to compare their station-keeping performances. For each simulation with the given nominal LPOs with the SRP, 100 sample points with different α , initial insertion errors and practical constraint parameters (see Table 2) are chosen. Table 6 lists the simulation results with λ=0.01, 0.0187, 0.03, and 0.04 m2 /kg. From this table, we find that although based on Eqs. (3) and (5) the SRP acceleration increases with the increase of λ, the average fuel cost  v of station keeping has no significant changes. We suspect that this is because we use the halo orbit with the corresponding λ as the corresponding nominal LPO. Therefore, the influence of the SRP on the stationkeeping problem is reduced. However, the average idle time T idle decreases with the increase of λ, which is a defect of our control strategy. Based on the data in Table 6, we conclude that our control strategy with a dead-band scheme can be applied to the stationkeeping problem with a large λ, and the nominal LPO with the SRP in the ephemeris model can efficiently reduce the influence of the SRP, especially for the fuel cost.

Monte-Carlo results

Halo

Vertical Lyapunov

Lissajous

Avg. v (m/s) Max. v (m/s) Min. v (m/s) Std. v (m/s) Avg.T idle (days) Max.T idle (days) Min.T idle (days) Std.T idle (days)

33.8848 40.0445 25.5184 3.8804 95.0288 133.7025 69.4111 18.5523

32.2608 35.4196 27.5736 2.4288 94.0736 130.7049 69.4111 15.7330

31.7329 41.8669 26.5326 3.9441 80.3524 109.2450 60.7977 11.2091

Simulation results for LPOs with longer durations are also studied. For example, the control strategy for the Lissajous orbit with the SRP in 2 years is illustrated in Fig. 16. Numerical calculation indicates that its  v is 36.1378 m/s and its T idle is 82.1553 days. By looking at the figure, we find the control strategy for 2 years is similar to that in one year: the control strategy of station keeping can track the nominal Lissajous orbit, and the position deviation to the nominal orbit dreal is maintained in the range of about 5∼30 km. The Monte-Carlo simulations using different nominal LPOs with the SRP in 2 years are listed in Table 7. In the Monte-Carlo simulations, for each simulation with the given nominal LPOs with the SRP, 100 sample points with different α , initial insertion errors and practical constraint parameters are chosen. Similar to the conclusion in Table 5, we also find that the fuel cost  v and the total idle time T idle are similar among different LPOs. Therefore, we infer that our control strategy could be applied to the long-term station keeping for different translunar LPOs. Station-keeping performances for different LPOs with the SRP in different simulation time are shown in Fig. 17, where points denote the Monte-Carlo results for different LPOs and solid lines link average values of corresponding Monte-Carlo simulations. In these figures, results of different LPOs are similar for the given simulation time. When the simulation time is larger than 0.5 years, both  v and T idle seem to grow linearly as the simulation time are added. This phenomenon is what one would intuitively expect. As we can see from in Fig. 16, when the time t is larger than 0.5 years, the station-keeping strategy is in steady-state, well beyond the initial transient period. For purely steady-state operation it would be expected that both  v and T idle increase linearly with simulation time. Therefore, we postulate that our control strategy under practical constraints can be applied to the long-term station keeping of translunar LPOs, and its performances,  v and T idle , have approximate linear relationships with the duration time of station keeping. Based on this conclusion, we infer that average  vs of the halo orbit, vertical Lyapunov orbit, and Lissajous orbit, for 10 years are about 156.1617 m/s, 158.5811 m/s, and 142.6511 m/s, respectively; average T idle s of halo orbit, vertical Lyapunov, and Lissajous orbit, for 10 years are about 493.8819 days, 502.0284 days, and 375.8150 days, respectively. Obviously, station keeping of the Lissajous orbit needs less fuel cost, but its idle time is also smaller. Using the impulsive maneuver and the DSMC, Lian et al.

14

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

Fig. 16. Control strategy for the Lissajous orbit in 2 years: (a) position deviation to the nominal orbit, (b) control input, (c) controlled Lissajous orbit in the Earth-Moon rotating frame, where red and green lines denote the idle and thrust segments, respectively.

Fig. 17. (a) Position distance of Dead-band control, (b) thrust of Dead-band control.

[9] investigated station keeping of translunar LPOs in 10 years. Based on their results, under the similar parameters of practical constraints in our simulations, the total fuel cost of station keeping for halo orbits are about 169 m/s, and the total fuel cost for Lissajous orbits are also smaller, about 150∼156 m/s. Hence, compared with the result of Lian et al. [9], our station-keeping strategy needs less fuel costs. Last but not the least, we pay attention to the influence of errors with different distributions on station keeping. For the sake of space, we only take into consideration different initial insertion errors, and other errors and practical constraints are kept as the data in Table 2. For the initial insertion error, we choose two kinds of distributions herein: the normal distribution and the uniform distribution. In addition, the SRP is also considered in the simulation, and we choose the halo orbit with the SRP in Fig. 2 as the nominal orbit. The simulation time is one year. In the Monte-Carlo simulations, for the given nominal halo orbits, 100 sample points with different α , initial insertion errors and practical constraint parameters are chosen. The simulation results of the initial insertion error with the normal distribution have been obtained in Table 5, so we only need to calculate the case of the uniform distribution. We assume that the two kinds of distributions of the initial insertion error have the same mean and the same standard deviation. Because the initial insertion errors are assumed to have standard normal distributions with zero means and standard deviations of σ X ,Y , Z = 100 km and σ X˙ ,Y˙ , Z˙ = 1 cm/s in previous simulations, so when the initial insertion errors have uniform distributions, their maximum and minimum initial position errors should be 173.2051 and −173.2051 km, respectively, and their maximum and minimum initial velocity errors should be 1.7321 and −1.7321 cm/s, respectively. Based on the above assumption, the simulation results of the case of the uniform distribution are listed in Table 8. As we can see from the table, although the initial insertion errors have different distributions, their average fuel costs and average idle time have no significant differences. Therefore, our control strategy can be applied to other station-keeping problems with different distributions.

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

Table 8 Comparison of station keeping with different initial insertion errors. Monte-Carlo results

Normal distribution

Uniform distribution

Avg. v (m/s) Max. v (m/s) Min. v (m/s) Std. v (m/s) Avg.T idle (days) Max.T idle (days) Min.T idle (days) Std.T idle (days)

18.6002 23.6432 14.8285 2.9872 45.1722 69.9073 24.8129 12.9095

18.5144 23.4524 15.2729 2.3338 41.5628 57.2297 30.1103 7.5151

15

Declaration of competing interest No competing interest. Acknowledgements This work is supported by the Natural Sciences and Engineering Research Council of Canada through a Discovery Accelerator Supplement under Grant RGPAS-493042-2016. References

5. Conclusions In this paper, the station-keeping strategy for translunar libration point orbits (LPOs) was investigated in the ephemeris model using the continuous thrust. Ephemeris models with and without the solar radiation pressure (SRP) acceleration were proposed for the numerical simulation. Three kinds of translunar LPOs, including halo orbits, vertical Lyapunov orbits, and Lissajous orbits, were constructed in ephemeris models with and without the SRP as nominal orbits of our research. Numerical results indicated between LPOs with and without the SRP, both the position difference and the velocity difference cannot be ignored in the stationkeeping problem. A backstepping technique with the dead-band scheme was extended to station keeping in the ephemeris model. The validity and rapidity of this control strategy were verified by the numerical computation. Furthermore, optimal parameters of the dead-band scheme were discussed by the Monte-Carlo simulations. Numerical results indicated that the dead-band scheme with an upper boundary of 20 km and a lower boundary of 5 km was a well balanced option between the fuel cost and the idle time. Station keeping under practical constraints, caused by the navigation system and the executive system, was studied in ephemeris models with and without the SRP. A control strategy with practical constraints was presented and applied to station keeping of LPOs. In numerical simulations without the SRP, we found that using just few more fuel costs, the control strategy with a dead-band scheme could have more than twice idle time than that without a deadband scheme. Hence, we postulated that the dead-band scheme was preferable in station keeping. In addition, the numerical calculation showed that the existence of navigation errors could result in less idle time but had no obvious influences on the fuel cost. Furthermore, the SRP was taken into account in station keeping for LPOs. Based on numerical simulations, we preferred LPOs obtained in the ephemeris model with the SRP to those without the SRP as nominal orbits, because the former had less fuel costs and more idle times. The investigation of the navigation interval time indicated that a shorter interval time could result in less fuel costs and more idle times, but in practical missions, this value was limited by the navigation system. The influence of the area-to-mass ratio of the spacecraft on station keeping was discussed based on the Monte-Carlo simulations. We found our control strategy with a dead-band scheme can be applied to the station-keeping problem with a large area-to-mass ratio, and the nominal LPO with the SRP in the ephemeris model can efficiently reduce the influence of the SRP, especially for the fuel cost. More Monte-Carlo simulations testified that our control strategy could be applied to the long-term station keeping for different translunar LPOs. The fuel cost and the idle time have approximate linear relationships with the duration time of station keeping. In addition, simulation results also showed that our control strategy can be applied to other station-keeping problems with different distributions. Results obtained in this paper provide useful suggestions for station keeping using the continuous thrust in future.

[1] ISECG, The Global Exploration Roadmap, Tech. rep., National Aeronautics and Space Administration, 2011, https://permanent.access.gpo.gov/gpo21917/ 591067main_GER_2011_small.pdf. [2] V. Angelopoulos, The ARTEMIS mission, in: The ARTEMIS Mission, Springer, 2010, pp. 3–25. [3] Y. Lian, G. Gómez, J.J. Masdemont, G. Tang, A note on the dynamics around the Lagrange collinear points of the earth–moon system in a complete solar system model, Celest. Mech. Dyn. Astron. 115 (2) (2013) 185–211, https://doi. org/10.1007/s10569-012-9459-2. [4] W. Wiesel, W. Shelton, Modal control of an unstable periodic orbit, J. Astronaut. Sci. 31 (1983) 63–76. [5] C. Simó, G. Gómez, J. Llibre, R. Martinez, J. Rodriguez, On the optimal station keeping control of halo orbits, Acta Astronaut. 15 (6–7) (1987) 391–397, https://doi.org/10.1016/0094-5765(87)90175-5. [6] G. Gomez, K. Howell, J. Masdemont, C. Simo, Station-keeping strategies for translunar libration point orbits, in: AAS/AIAA Spaceflight Mechanics Conference, Monterey, CA, 1998, paper AAS 98-168. [7] X. Hou, L. Liu, J. Tang, Station-keeping of small amplitude motions around the collinear libration point in the real Earth–Moon system, Adv. Space Res. 47 (7) (2011) 1127–1134, https://doi.org/10.1016/j.asr.2010.12.005. [8] K. Howell, H. Pernicka, Station-keeping method for libration point trajectories, J. Guid. Control Dyn. 16 (1) (1993) 151–159, https://doi.org/10.2514/3.11440. [9] Y. Lian, G. Gomez, J.J. Masdemont, G. Tang, Station-keeping of real Earth–Moon libration point orbits using discrete-time sliding mode control, Commun. Nonlinear Sci. Numer. Simul. 19 (10) (2014) 3792–3807, https://doi.org/10.1016/j. cnsns.2014.03.026. [10] X. Bai, J.L. Junkins, Modified Chebyshev-Picard iteration methods for stationkeeping of translunar halo orbits, Math. Probl. Eng. (2012) 926158, https:// doi.org/10.1155/2012/926158. [11] M. Shirobokov, S. Trofimov, M. Ovchinnikov, Survey of station-keeping techniques for libration point orbits, J. Guid. Control Dyn. 40 (5) (2017) 1085–1105, https://doi.org/10.2514/1.G001850. [12] J.V. Breakwell, A.A. Kamel, M.J. Ratner, Station-keeping for a translunar communication station, Celest. Mech. 10 (3) (1974) 357–373, https://doi.org/10.1007/ BF01586864. [13] S. Jia, Y. Jia, S. Xu, Q. Hu, Maneuver and active vibration suppression of freeflying space robot, IEEE Trans. Aerosp. Electron. Syst. 54 (3) (2017) 1115–1134, https://doi.org/10.1109/TAES.2017.2775780. [14] S. Jia, J. Shan, Optimal actuator placement for constrained gyroelastic beam considering control spillover, J. Guid. Control Dyn. (2018) 1–9, https://doi.org/ 10.2514/1.G003560. [15] V. Deshmukh, S. Sinha, Control of dynamic systems with time-periodic coefficients via the Lyapunov-Floquet transformation and backstepping technique, Modal Anal. 10 (10) (2004) 1517–1533, https://doi.org/10.1177/ 1077546304042064. [16] M. Nazari, W.M. Anthony, E. Butcher, Continuous thrust stationkeeping in Earth-Moon L1 halo orbits based on LQR control and Floquet theory, in: AIAA/AAS Astrodynamics Specialist Conference, San Diego, California, 2014, p. 4140. [17] Y. Ulybyshev, Long-term station keeping of space station in lunar halo orbits, J. Guid. Control Dyn. 38 (6) (2014) 1063–1070, https://doi.org/10.2514/ 1.G000242. [18] A. Narula, J.D. Biggs, Fault-tolerant station-keeping on libration point orbits, J. Guid. Control Dyn. (2017) 1–9, https://doi.org/10.2514/1.G003115. [19] H. Peng, Y. Liao, X. Bai, S. Xu, Maintenance of libration point orbit in elliptic Sun–Mercury model, IEEE Trans. Aerosp. Electron. Syst. 54 (1) (2018) 144–158, https://doi.org/10.1109/TAES.2017.2739938. [20] W.M. Folkner, J.G. Williams, D.H. Boggs, R.S. Park, P. Kuchynka, The Planetary and Lunar Ephemerides DE430 and DE431, Interplanetary Network Progress Report (196), 2014, pp. 1–81. [21] J. Bookless, C. McInnes, Control of Lagrange point orbits using solar sail propulsion, Acta Astronaut. 62 (2–3) (2008) 159–176, https://doi.org/10.1016/j. actaastro.2006.12.051. [22] V. Domingo, B. Fleck, A.I. Poland, The SOHO mission: an overview, Sol. Phys. 162 (1–2) (1995) 1–37, https://doi.org/10.1007/BF00733425.

16

Y. Qi, A. de Ruiter / Aerospace Science and Technology 94 (2019) 105376

[23] K.C. Howell, Three-dimensional, periodic, ‘halo’ orbits, Celest. Mech. 32 (1) (1984) 53–71, https://doi.org/10.1007/BF01358403. [24] J.S. Parker, R.L. Anderson, Low-Energy Lunar Trajectory Design, John Wiley & Sons, 2014. [25] X. Hou, L. Liu, On quasi-periodic motions around the collinear libration points in the real Earth–Moon system, Celest. Mech. Dyn. Astron. 110 (1) (2011) 71, https://doi.org/10.1007/s10569-011-9340-8. [26] D.L. Richardson, Analytic construction of periodic orbits about the collinear points, Celest. Mech. 22 (3) (1980) 241–253, https://doi.org/10.1007/ BF01229511.

[27] B.G. Marchand, K.C. Howell, R.S. Wilson, Improved corrections process for constrained trajectory design in the n-body problem, J. Spacecr. Rockets 44 (4) (2007) 884–897, https://doi.org/10.2514/1.27205. [28] T.H. Sweetser, S.B. Broschart, V. Angelopoulos, G.J. Whiffen, D.C. Folta, M.-K. Chung, S.J. Hatch, M.A. Woodard, Artemis mission design, Space Sci. Rev. 165 (1–4) (2011) 27–57, https://doi.org/10.1007/s11214-012-9869-1. [29] D. Cielaszyk, B. Wie, New approach to halo orbit determination and control, J. Guid. Control Dyn. 19 (2) (1996) 266–273, https://doi.org/10.2514/3.21614.