- Email: [email protected]
Phasing Delta-V for transfers from Sun–Earth halo orbits to the Moon
Acta Astronautica 127 (2016) 464–473
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/aa
Phasing Delta-V for transfers from Sun–Earth halo orbits to the Moon Hongru Chen a,n, Yasuhiro Kawakatsu b, Toshiya Hanada a a b
Department of Aeronautics and Astronautics, Kyushu University, Fukuoka, Japan Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency (ISAS/JAXA), Sagamihara, Japan
art ic l e i nf o
a b s t r a c t
Article history: Received 8 June 2015 Received in revised form 9 December 2015 Accepted 3 May 2016 Available online 26 May 2016
Inspired by successful extended missions such as the ISEE-3, an investigation for the extended mission that involves a lunar encounter following a Sun-Earth halo orbit mission is considered valuable. Most previous studies present the orbit-to-orbit transfers where the lunar phase is not considered. Intended for extended missions, the present work aims to solve for the minimum phasing ΔV for various initial lunar phases. Due to the solution multiplicity of the two-point boundary value problem, the general constrained optimization algorithm that does not identify multiple feasible solutions is shown to miss minima. A two-step differential corrector with a two-body Lambert solver is developed for identifying multiple solutions. The minimum ΔV associated with the short-way and long-way approaches can be recovered. It is acquired that the required ΔV to cover all initial lunar phases is around 45 m/s for the halo orbit with out-of-plane amplitude Az greater than 3.5×105 km, and 14 m/s for a small halo orbit with Az ¼ 1×105 km. In addition, the paper discusses the phasing planning based on the ΔV result and the shift of lunar phase with halo orbit revolution. & 2016 IAA. Published by Elsevier Ltd. All rights reserved.
Keywords: Extended missions Circular restricted three-body problem Halo orbits Lunar transfers Two-point boundary value problem
1. Introduction The halo orbit about a libration point in a three-body system has features such as relatively constant distances and orientation with respect to the two primary bodies, which are favorable for scientific observation, thermal control and communication. A famous mission is the ISEE-3. ISEE-3 achieved several missions including monitoring solar wind in a Sun–Earth halo orbit, exploring geo-tail and visiting comets beyond original expectations [1]. The success of the mission is owing to the elaborate planning of transfers where the three-body dynamics and multiple lunar gravity assists were used. There have been several libration point missions in addition to the ISEE-3, and many are being planned. The success of ISEE-3 and the increasing interest in halo orbit missions give rise to the need of investigation of extended missions following them. Utilizing the low-energy transfer along the manifold associated with a libration point orbit, various missions can be developed. One option for the extended mission of a Sun–Earth halo orbit mission is going to the Moon. Following the success of the HITEN mission, several studies and mission design methods have demonstrated the low-energy transfers to the Moon under the gravitational influences of the Sun and the Moon [2–6]. The extended mission of THEMIS, ARTEMIS, reached to the Sun–Earth L1 n
Corresponding author. E-mail address: [email protected] (H. Chen).
http://dx.doi.org/10.1016/j.actaastro.2016.05.003 0094-5765/& 2016 IAA. Published by Elsevier Ltd. All rights reserved.
region and then an Earth-Moon L2 Lissajous orbit applying manifold dynamics [7]. Another promising option is the interplanetary exploration. However, low-energy escape only permits accessing a few asteroids in a narrow zone. For instance, the low-energy transfer from the Earth to the Mars needs an additional ΔV of around 2 km/s [8,9]. On the other hand, lunar gravity assists have proven to be an effective way to increase escape energy [10–14]. The ISEE-3 transferred from the Sun–Earth L1 halo orbit to the Moon and then performed multiple lunar gravity assists for Earth escape. A previous study investigates the combination of the unstable manifold and multiple lunar swingbys for Earth escape [15]. The key remarks obtained from that study are: 1) manifold-guided lunar swingbys can achieve much higher characteristic energy (C3) with respect to the Earth than direct escape along the manifold; 2) for the case that the V1 with respect to the Moon is not sufficient for high-energy Earth escape, a second lunar swingby can efficiently increase the C3 to a level of 2.7 km2/s2 in 70–90 days. Transfers to the Moon utilizing manifold dynamics have been intensively studied and well demonstrated [2–6]. However, previous studies are to find low-cost transfer paths and meanwhile decide the optimal departure orbit phase or terminal lunar phase. Using numerical techniques, Ref. [4] solves the zero-ΔV transfer from a Sun–Earth L2 halo orbit to an Earth–Moon L1 halo orbit. Different from those studies on the orbit-to-orbit transfers, the motivation of this paper is the extended mission where the departure orbit, i.e. a Sun–Earth halo orbit, is not pre-phased for a lunar encounter. The previous work (Ref. [15]) assumes the transfers to the Moon are also ballistic. The present paper is a
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
follow-up work concerning a phasing maneuver for encountering the Moon. Some phasing ΔV were used to target the Moon in ISEE3 and ARTEMIS. Obtaining the required phasing ΔV for different conditions will facilitate the design of similar extended missions that involve a lunar encounter. The study is also motivated by JAXA's mission DESTINY (Demonstration and Experiment of Space Technology for INterplanetary voYage) [16]. DESTINY will be launched to a halo orbit about the Sun–Earth L2 point in 2021, and the possibility of extending the mission to visit a heliocentric body is currently under discussion. Section 2 presents the basis of the circular restricted threebody problem. Section 3 gives the overview and parameterization of the phasing problem for the transfer from a Sun–Earth halo orbit to the Moon. Section 4 presents the computation procedures and the results of minimum required phasing ΔV for various initial lunar phases. In attempts of optimization, the problem of missing minima is found due to the solution multiplicity of the two-point boundary value problem (TPBVP). To tackle this problem, a sequential differential corrector aided by a two-body Lambert solver is developed to identify the multiple solutions to a TPBVP robustly. The required ΔV to cover all initial lunar phases is acquired. In addition, the characteristics of the transfer trajectories with minor phasing ΔV are discussed for different mission scenarios, namely, lunar mission and lunar gravity-assisted Earth escape. Section 5 discusses the phasing planning based on the result of minimum ΔV, and the relationship between the initial lunar phase and the halo orbit revolution. Conclusions are given in the final section.
2. System dynamics 2.1. Circular restricted three-body problem Trajectories are considered in the circular restricted three-body problem (CR3BP) of the Sun–Earth system. The CR3BP assumes two primary bodies m1 and m2 moving in a circular orbit about their barycenter. The mass of the third body is negligible compared to the masses of the two primaries. The two primaries can be the Sun and Earth, the Earth and Moon or etc. The rotating coordinate system with the barycenter at the origin and the primaries fixed on the x-axis is chosen to describe the motion of the third body. For convenience, the angular velocity of the rotating frame, the total mass and the distance between the two primaries are normalized to 1. μ is the ratio of the mass of m2 to the mass of m1. Then, masses of m1 and m2 become 1-μ and μ in the normalized model. The coordinates of m1 and m2 are [ μ, 0, 0] and [1 μ, 0, 0], respectively, as shown in Fig. 1. The equations of motion of the third body are:
x¨ − 2ẏ = ∂U /∂x
(1)
y¨ + 2ẋ = ∂U /∂y
(2)
z¨ = ∂U /∂z
(3)
where the gravitational potential U in the system is
U = (x2 + y2 )/2 + (1 − μ)/rP1 + μ /r P2
(4)
where
r P1 =
(x + μ)2 + y2 + z2 , r P2 =
(x − 1 + μ)2 + y2 + z2
(5)
465
Fig. 1. Circular restricted three-body system and Lagrangian points.
from the equations of motion. They are also referred to as the libration points or Lagrangian points labeled as L1,…, L5. The three collinear points L1, L2 and L3 are unstable. The geometries of the five equilibrium points are depicted in Fig. 1. 2.2. Halo orbits and invariant manifolds The linearized equations at the collinear equilibrium points reveal that the motion in the vicinity of these points consists of the in-plane and out-of-plane oscillatory modes. A halo orbit is a periodic orbit with a sufficiently large out-of-plane amplitude Az such that the out-of-plane oscillatory frequency becomes equal to the in-plane oscillatory frequency [18]. Halo orbits are generally computed numerically [19]. The monodromy matrix is the state transition matrix for one period of a periodic orbit. By examining the eigenvalues and eigenvectors of the monodromy matrix, the insight of the orbit dynamics can be revealed. For halo orbits, there is a pair of eigenvalues with λ1 λ2 ¼ 1 [18]. The dominant eigenvalue λ1 is around 1500 for Sun–Earth L1/L2 halo orbits. A small displacement along the eigenvector of the dominant eigenvalue (divergent eigenvector for short) will propagate into exponentially increasing divergence from the halo orbit quickly, since the displacement will grow λ1 times in a period. Such a departure trajectory follows the unstable invariant manifold. A derived manifold trajectory can go earthward or anti-earthward. On the other hand, impulsive changes of both position and velocity are not practical in real missions. Δv are added to generate manifold trajectories instead. The Δv direction corresponding to the divergent eigenvector is computed based on Appendix A. To globalize the manifold, small Δv in the divergent directions are added to a series of points along the halo orbit. The propagation of the perturbed states yields an unstable manifold with a tube structure, as is shown in Fig. 2. Conversely, backward propagating the states perturbed by small Δv in convergent directions corresponding to the eigenvalue smaller than 1 (λ2) will result in the stable manifold that asymptotically converges into the halo orbit.
3. Phasing problem
An energy quantity,
J = 2U − (x2̇ + y2̇ + z2̇ )
(6)
holds in the system. It is known as the Jacobi integral. Ref. [17] presents the detail of resolving the equilibrium points
As both of the unstable invariant manifold and lunar gravity assist are useful in delivering a spacecraft to deep space, our previous study applies the combination of them for Earth escape. That study shows this strategy can achieve a considerable C3 for
466
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Az = 4x105 km
−3
6
x 10
corrected trajectory ωΔtf
Δxf, Δtf
rM (tf) xf, tf
4 1
initial trajectory (along manifold)
y ,AU
2 0
2
L2
4
Moon’s orbit
3
x 1 -, t 1 -
x1+, t1+ Δv
−2 Moon’s orbit
−4 −6
Fig. 3. Correct a manifold trajectory to target the Moon.
0.995
1
1.005
1.01
x ,AU −3
x 10
z, AU
2 2
0
3
1
L2
4
−2 0.995
1
1.005
1.01
x, AU Fig. 2. The intersection of the manifold tube associated with the halo orbit with Az ¼ 4×105 km and Moon's orbit exhibits four types of intersection.
interplanetary missions. An overview of the previous study is presented in this section, which leads to the discussion on the phasing problem. 3.1. Background The spacecraft is assumed initially in a Sun–Earth L2 halo orbit. The inclination and eccentricity of the Moon's orbit are small and assumed zero. This simplified model is sufficient for acquiring reliable information for practical mission design [20]. Fig. 2 shows the earthward trajectories of the unstable manifold of the halo orbits with Az ¼4×105 km and the Moon's orbit. The grey circles mark the intersections of the manifold trajectories and the ecliptic plane, outlining the intersection lines of the manifold tube and the ecliptic plane. It is found that the manifold and the Moon's orbit intersect at four points after the first crossing through Moon’s orbit. These four types of intersection exist for the wide Az range investigated (from 1×105 km to 5×105 km). The four types of intersection are numbered in order of the lunar phases. The previous study shows, if the Moon is exactly at the intersection points at the manifold arrival, a following single or double lunar swingby can lift the C3 energy with respect to the Earth to around 2.7 km2/s2, which is far greater than a direct escape from the halo orbit via the anti-earthward unstable manifold. However, since the manifold tube is fixed for a specified halo orbit, the positions of the four intersection points are fixed. The Moon may not be at the intersection points at the time the manifold reaches. Therefore, a phasing maneuver to achieve a lunar encounter is considered, and the cost of it is investigated in this paper. 3.2. Parameterization Define the epoch the spacecraft passes a reference point rref in the halo orbit as the reference epoch tref. rref is set at [x L2 + Ax , 0, Az ] in this study, where [x L2, 0, 0] is the coordinate of L2
and Ax is the x-amplitude of the orbit. Define the lunar phase at tref as the initial lunar phase θ0. θ0 is related with the launch date and halo orbit revolution. The halo orbit size Az, which is associated with the Jacobi integral J is known once the mission is decided. For a given Az, the required ΔV to target the Moon for various θ0 would be a useful reference for evaluating and planning the extended mission involving a lunar encounter. It is reasonable to assume that the scheme of leaving the halo orbit via the manifold and then applying a mid-course ΔV to target the Moon is an optimal choice. To find the optimal patch point along the manifold, the manifold is parameterized. Two time variables τ and t1 are introduced. τ is the time of the manifold insertion point from tref. The magnitude of manifold insertion Δv is set to be 1 cm/s, then, τ determines a manifold trajectory. t1 is the propagation time of the manifold and also the timing of the phasing maneuver. The final epoch tf determines the lunar encounter position rM, which is a function of θ0 and τ þtf. The maneuver at t1 corrects the manifold trajectory to reach rM at tf. Fig. 3 schematically depicts the phasing problem. A control of V1 with respect to the Moon can further enhance the mission design. The V1 is a critical factor to the mission using lunar gravity assists. However, for a lunar mission, the V1 should be reduced. The control of V1 should depend on the specific mission scenario. Discussions for various mission scenarios are beyond the scope of this paper. As for using lunar gravity assists to achieve Earth escape, studies have shown that even if the V1 is small, a second lunar swingby without additional ΔV cost can efficiently lift the escape energy to the theoretical maximum level in 2-3 months [14]. Therefore, the objective of phasing in this study is simply to achieve a lunar encounter with minimum ΔV. However, characteristics of the lunar transfer trajectories with minor phasing ΔV will be revealed and discussed for different mission scenarios in the subsequent section.
4. Optimization of phasing ΔV 4.1. Constrained optimization Trajectories are subject to Eqs. (1)–(3). For given Az and optimization is to find the 6-element set, T y = ⎡⎣ τ , t1, tf , ΔvT⎤⎦
θ0, the (7)
that satisfies the equality boundary constraint,
c = r M (tf ) − rf = 0
(8)
where rf is the position of the spacecraft at tf, and at the same time minimizes the cost function,
ΔV = Δv
(9)
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
467
A number of initial guesses of y should be tried in order to find the global minimum. It is prescribed that,
−T < τ < 0
(10)
where T is the period of the halo orbit, which is a function of Az. This prescription suggests a slot of one halo orbit period before tref for targeting a lunar encounter. 20 τ′ are sampled evenly within the prescribed slot as the initial guesses. (Prime and asterisk superscripts are used to indicate initial guesses and optimal solutions, respectively.) Then, the full range of manifold trajectories is taken into account. Note that, once τ* is solved, it does not mean the spacecraft must be inserted to the manifold at τ*. If the spacecraft misses τ*, it can be inserted to the same manifold trajectory at a later τ by paying a correspondingly larger manifold insertion ΔV. As mentioned in Section 3.1, encounter chances and closest approaches can occur at every intersection of the manifold tube and the Moon's orbit. There can be infinite local minima if tf is not constrained. To simplify the problem, this study only concerns the first encounter chances. Then, tf′ is chosen to be the epoch of the first closest approach to the Moon along the manifold trajectory. As for t1′, through experiments with some boundary conditions, it is observed that there are usually two local minima of ΔV along the t1 coordinate (an example will be presented in Fig. 7). t1′ are chosen to be 30 and 120 days prior to tf′ based on the observation. It becomes difficult for the numerical computation to converge as t1′ approaches 0, because of the high sensitivity in the divergent direction in the halo orbit regime (λ1 E 1500 compared to other eigenvalues r 1). Therefore, the ranges and sampled initial guesses of τ, t1 and tf are considered reasonable and sufficient for the problem setting. Δv′ is 0. The MATLAB constrained optimization solver, fmincon, is used in this study. Minimum ΔV is computed for the θ0 varied from 0° to 360°. θ0 is sampled every 5°. Eventually, the required ΔV budget to cover the full range of θ0 can be obtained, providing a reference for extended missions. Fig. 4 summarizes the procedures to generate the minimum phasing ΔV as a function of θ0. 4.2. Results and discussions of constrained optimization The halo orbit with Az ¼4×105 km is taken as an example to investigate the minimum phasing ΔV, ΔVmin, as a function of θ0.
Fig. 5. The minimum phasing ΔV and the optimal τ solved by the constrained optimization (Az ¼ 4×105 km).
The ΔVmin vs θ0 computed by the workflow in Fig. 4 is shown in Fig. 5. The ΔVmin profile is shown made up of several profiles. Let τ1,…, τ4 represent the τ of the four manifold trajectories that intersect with the Moon's orbit, and θ01…, θ04 correspondingly represent the θ0 conditions under which the four manifold trajectories can meet the Moon without phasing ΔV. The τ* of each profile are found near τ1,…, τ4 for the θ0 near θ01…, θ04 (see the τ* variation measured by the scale to the right), respectively. Then, the minima can be grouped into four groups corresponding to the four types of intersection. To get an insight, Fig. 6 shows the ΔVmin variations associated with Type 1 and Type 2 in the θ0 interval [300° 340°]. It can be seen that the ΔVmin associated with the two types take turns to stand out (i.e., to have lower ΔV for the θ0 by comparison) in this interval. Because a sufficient number of τ′ are tried in the computation and the τ of the two types are quite different, the program can successfully identify the minima of the two groups and output the lower ones. On the other hand, it is observed that the zero minimum of the Type-3 intersection as well as the minima associated with it are missing in the result. This problem is found related with the solution multiplicity of the two-point boundary value problem (TPBVP). It is well known that, in the two-body problem, there are multiple solutions to the TPBVP, namely, the short-way, long-way and multiple-revolution approaches. The manifold trajectory of the Type-3 intersection can be regarded as a long-way approach in the near-Earth realm; while the other three are short-way approaches (see Fig. 2). For the group-3 and group-4 manifold trajectories, the required phasing ΔV that leads to short-way and long-way approaches may not differ greatly. For illustration, Fig. 7 shows the ΔV variations associated with the two approaches for various two points, which are defined by sets of t1 and tf, for a given θ0 near the θ03 and a given τ near τ3. It can be seen that, each solution profile exhibits two local minima. Moreover, for a considerable range, the ΔV for the two approaches are close to each other. The constrained optimization algorithm solves TPBVP implicitly. Once the optimizer converges to a feasible solution, it is inevitable for the gradient-based approach to be trapped in the current solution regime. Therefore, the lowest ΔV minimum can be possibly missed by a general constrained optimizer. One may try multiple initial guesses for Δv. However, as there are actually three variables in Δv, if no good knowledge on their ranges is available, it may take excessively long time on many random initial guesses in order to guarantee all feasible solutions. Good estimates of the Δv of both trajectory motions are desirable in this problem. A TPBVP solver is developed in the next subsection for the purpose of distinguishing different solutions. 4.3. Sequential differential corrector with a Lambert solver for TPBVP
Fig. 4. Workflow of computing the minimum phasing ΔV vs θ0 for a given Az using constrained optimization.
To identify the Δv for the desired trajectory motion, the first step is to obtain a close guess utilizing a two-body Lambert solver
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Group 1 Group 2 (near τ2 = -144 days) (near τ1 = -17.7 days)
468
ΔVmin(group 1) > ΔVmin(group 2)
ΔVmin(group 1) < ΔVmin(group 2)
ΔVmin(group 1) > ΔVmin(group 2)
Fig. 6. The minimum ΔV of group 1 and group 2 for the θ0 varied from 300° to 340° (Blank areas are where the ΔV is too high to be displayed or failed to be solved by the program).
Fig. 7. ΔV contours of the short-way and long-way approaches. (Blank areas are where ΔV is too high to be displayed or failed to be solved by the program).
that can identify multiple solutions in the two-body problem. The trajectory is segmented at an auxiliary epoch t2 when the manifold trajectory enters a defined Earth-centered sphere with a radius of 8×105 km. The sphere is considered as the boundary between the three-body and two-body realms. Then, the t1 t2 arc is seen as a three-body trajectory, while the t2 tf arc a two-body trajectory. The t1 t2 arc is propagated, which gives the position r2 and velocity v2 at t2. Let "i" superscript indicate vectors referenced to the inertial frame. The t2 tf arc connecting r2 and rfi is solved by a two-body Lambert solver. rfi is the rf rotated by an angle of Ω (tf t2), where Ω is the rotating rate of the Sun Earth system. The Lambert solver used in this study is based on Bate's method [21] and can identify the solution (v2i+) for a desired (short-way or long-way) trajectory motion. Then, the corresponding v2+ referenced to rotating frame can be obtained. There appears an inconsistency,
Fig. 8. Break the trajectory at an auxiliary epoch t2, and then smoothen the connection of the two-body and three-body arcs by differential correction.
Δv2 = v 2+ − v 2−
(11)
By reducing Δv2, a good initial guess for the Δv to establish the desired trajectory can be obtained. The v2 can be tuned by a differential correction based on the propagated state transition matrix of the t1 t2 arc. The method of deriving the state transition matrix for CR3BP can be found in Ref. [19]. A recent study [22] gives the analytic expressions of the partial derivatives of the twobody Lambert problem, providing the gradients to correct v2 þ as well. Fig. 8 illustrates the correction concept. In detail, for the t1 t2 arc, differentiating the state x2 at t2 with respect to the state x1 at t1 yields,
δ x2 = Φ (t2, t1) δ x1
(12)
where Φ(t2, t1) is the 6×6 state transition matrix that maps the state at t1 to the state at t2. Eq. (12) can be expanded to
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
⎡ δr2 ⎤ ⎡ Φ rr (3 × 3) Φ rv (3 × 3) ⎤ ⎡ δr1 ⎤ ⎥⎢ ⎥ ⎢ ⎥=⎢ ⎣ δ v2 ⎦ ⎣ Φ vr (3 × 3) Φ vv (3 × 3) ⎦ ⎣ δ v1⎦
−3
x 10
(13)
5
((t2, t1) is dropped for brevity). Then, the partial derivative of v2 with respect to Δv is
4
∂v −2 = Φ vv ∂Δv
2
(15)
The r2 change results in two folds of v2 þ changes. One is the change of inertial velocity V 2i + computed by the Lambert solver. The other is the change of the frame-induced item Δv i2→ r ( r2), which is added to the inertial velocity to acquire the velocity referenced to the rotating frame; that is,
v 2+ = v i2+ + Δv 2i→ r
(16)
Δv i2→ r ( r2) is calculated from
Δv 2i→ r
⎡ 0 1 0⎤ ⎢ ⎥ = ⎢ − 1 0 0⎥ r2 = Rr2 ⎣0 0 0⎦
(17)
(18)
The complete derivation and expression for L are given in Ref. [22]. With Eqs. (15), (17) and (18), the derivative of Eq. (16) with respect to Δv can be written as,
∂v 2+ = (L + R) Φ rv ∂Δv
(19)
Then, integrating Eqs. (11), (14) and (19) results in the correction for Δv at the i-th iteration, which is
Δv (i) − Δv (i − 1) = (Φ vv − (L + R) Φ rv )−1Δv 2(i − 1)
y, AU
1 tf Moon
0 −1
short way
long way
t1
−2 −3 −4 −5 0.998
1
1.002
1.004
1.006
1.008
1.01
x, AU
Fig. 9. An example of the short-way and long-way approaches to the Moon identified by the TPBVP solver.
4.4. Unconstrained optimization with the developed TPBVP solver
A partial derivative of the Lambert problem is,
L = ∂v i2+/∂r2
3
(14)
The change of Δv will result in a change of r2, which is expressed as,
∂r2 = Φ rv ∂Δv
469
With the TPBVP solver developed to identify different trajectory motions, the boundary condition is always satisfied. Hence, the unconstrained optimization can be applied. The MATLAB unconstrained optimizer, fminunc, is used. As Δv is solved by the TPBVP solver, the number of the independent variables becomes three, which are, T y = ⎡⎣ τ , t1, tf ⎤⎦
(22)
Fig. 10 shows the procedures of using the developed TPBVP solver and unconstrained optimization. An effort was made to derive the first-order derivatives (G¼ ∂Δv/∂y) and second-order derivatives (H ¼∂2Δv/∂y2) to speed up the computation and ensure accuracy. The expressions of the first- and second-order derivatives are
(20)
Once Δv2 meets the tolerance (30 m/s in this study), Δv is considered a close estimate of the true solution. Note that there are situations that Δv2 cannot be reduced to the tolerance in many (100 is used in this work) iterations. A failure to converge suggests a large discontinuity between v2 and v2 þ , as well as a large ΔV for the desired trajectory. As the main interest is the minimum ΔV for given Az and θ0, for this case, the first step will assign an infinite value to the ΔV for the current guess, without going to the next step of correction. In addition, in the process of optimization, t1 may be tuned to be smaller than t2. In that case, a good Δv estimate is directly obtained from the two-body Lambert solver. The second step is correcting the Δv estimate obtained from the previous step, to acquire the accurate Δv for the three-body t1 tf trajectory. The trajectory with the estimated Δv is propagated to tf, and differential correction is carried out iteratively to gradually zero the position difference at tf (Eq. (8)) based on Φ(tf, t1); that is,
Δv (i) − Δv (i − 1) = − Φ rv (t f , t1)−1c
(21)
With the Δv guess obtained from the first step aided by the Lambert solver, the differential correction on the three-body trajectory will robustly converge to the Δv for the desired trajectory motion. Fig. 9 shows an example of the short-way and long-way approaches identified for a set of two points.
Fig. 10. Workflow of computing the minimum phasing ΔV vs θ0 for a given Az using unconstrained optimization with the TPBVP solver identifying short-way or longway solution.
470
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Fig. 11. The minimum phasing ΔV vs θ0 and the optimal τ solved by the three optimization programs (Az ¼4×105 km).
presented in Ref. [23]. On the other hand, it is time-consuming to tune the Δv starting from 0 for every y determined by the optimizer. After the Δv for a given set of two points, which are defined by y(i 1), is solved, the initial guess of the Δv for the next neighboring y(i) (if ||y(i) y(i 1)|| o1 day) is given by,
Δv (i) ′ = G (i − 1) (y (i) − y (i − 1) ) + Δv (i − 1)
(23)
If y(i) is not regarded as a neighboring point to the previous one, Δv′ starts from 0. In addition, to avoid repeated time-consuming integration to τ and t1, the manifold states are approximated by fitting the surface to the states on the grid of sampled τ and t1 using the bicubic convolution interpolation [24,25].The convolution interpolation function can also give derivatives with respect to τ and t1, with which G and H can be computed [23]. 4.5. Results and discussions of unconstrained optimization Fig. 11 displays the ΔVmin solved by the constrained optimization, unconstrained optimization with short-way and long-way solution controls, and the corresponding τ*. Optimizing with the trajectory motion control recovers the two corresponding ΔVmin variations. The results of constrained optimization generally follow the short-way solutions. The optimal short-way solutions are shown associated with the intersections of Type 1, Type 2 and Type 4, while the optimal long-way solutions are associated with the intersection of Type 3. Integrating the two ΔV profiles, it can be seen that a ΔV budget of 47 m/s can cover all initial lunar phases for a lunar encounter for Az ¼ 4×105 km. Fig. 12 shows the required ΔV to cover all θ0 as a
function of Az. The required ΔV decreases as Az decreases. In particular, the required ΔV decreases relatively quickly with Az starting from Az ¼ 3.5×105 km. By inspection, the profile of group 1 moves leftward, and that of group 2 moves rightward as Az decreases, while the profiles of group 3 and group 4 do not vary much. Group 1 starts to stand out for the most difficult θ0 instead of group 2 at Az ¼3.5×105 km. As ΔVmin peaks occur at the intersections of the ΔVmin profiles, the leftward movement of the profile of group 1 leads to a faster decreasing rate of the peak ΔVmin. For a small halo orbit with Az ¼1×105 km, a ΔV of 14 m/s is sufficient for targeting the Moon. The order of the ΔV cost is considered acceptable. For instance, ISEE-3 processed a ΔV capacity of 310 m/s after completing the original mission, and a ΔV of 72 m/s was allocated to carry out the transfer from the halo orbit to the escape trajectory [1]. Ref. [15] has presented the properties of the four types of ballistic lunar encounter for Earth escape. The V1 of the intersections of Type 3 and Type 4 with respect to the Moon is around 1.35 km2/ s2, which is much higher than that of Type 1 and Type 2. The corresponding lunar encounters of Type 3 and Type 4 can lead to effective lunar gravity-assisted Earth escape. Although the V1 of Type 1 and Type 2 are only around 0.4 km2/s2, a second lunar swingby can efficiently lift up the escape energy with some extra flight time. In addition, the manifold trajectories of Type 1 and Type 2 resemble the low-energy transfers to the Moon presented in Refs. [2–4]. Therefore, for an extended lunar mission there are two applicable options. One can solve for the ΔVmin vs θ0 associated with the Type-1 and Type-2 intersections, by optimizing near τ1 and τ2. Fig. 13 shows the results. Fig. 14 shows the corresponding phasing trajectories for various θ0. As the period of the halo orbit is roughly six times the lunar period, for a certain type of intersection, a corresponding fraction of manifold trajectories can be targeted to the Moon with minor ΔV (o100 m/s). A minor ΔV does not change the shape of a trajectory and the final lunar encounter situation too much, but mainly alter the time of flight for waiting for or catching up with the Moon. Fig. 13 also displays the corresponding V1 with respect to the Moon and ΔVmin þV1, which can serve as an estimate of the cost of a lunar mission. If only phasing ΔV is optimized, with the two options of transferring to the Moon, ΔVmin þ V1 varies from 475 m/s to 550 m/s for different θ0. The Moon's gravity might contribute to reducing the required ΔV in a more realistic model. However, computing the ΔVmin for Type-1 and Type-2 intersections can give a quick estimate of the required ΔV in preliminary lunar mission design. Fig. 15 displays the variations of τ*, t1* and tf*. τ*, t1* and tf* evolve
Fig. 12. The required ΔV to cover all θ0 as a function of Az.
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
x 10
471
Type 1
-3
6 manifold departure
4 phasing trajectories
Moon y, AU
2
L2
0 2 Δv
4 6 1
x 10
1.01
Type 2
-3
6
Fig. 13. The ΔVmin vs θ0 associated with the Type-1 and Type-2 lunar intersections, the resulting V1 with respect to the Moon and ΔVmin þ V1 (Az ¼ 4×105 km).
phasing trajectories
4
Δv
2 y, AU
linearly with θ0. The above characteristics result in the linear change of ΔVmin with θ0, as is shown in Figs. 11 and 13.
1.005 x, AU
L2
0 Moon
5. Phasing planning 2
If the manifold insertion slot in practical missions is planned in previous or next halo revolutions, the obtained ΔV vs θ0 is still (i ) refer to the epoch of rref of the i-th halo orbit reuseful. Let tref volution. The period of the halo orbit (Az ¼4×105 km) is 179.7 days and the period of the Moon's orbit (with respect to the Sun–Earth rotating frame) is 29.5 days. It can be inferred that the θ0 shifts 33° (i ) at every tref . Therefore, one can relate θ0 to different tref , and find the halo orbit revolution and corresponding θ0 that requires acceptably low phasing ΔV. For instance, supposing that a two-revolution slot is permitted for inserting the spacecraft to the (i − 1) (i ) , tref and manifold, the lunar phases at the consecutive epochs tref
4
manifold departure
6 1
1.005
1.01
x, AU Fig. 14. The transfer trajectories with minimum phasing ΔV for the Type-1 and Type-2 intersections. Arrows represent phasing Δv.
(i + 1) are θ0(i − 1) , θ0(i) and θ0(i + 1) , respectively. It should be noted that tref
the optimal manifold trajectory solved for θ0(i − 1) can be targeted (i − 1) (i ) during the revolution between tref and tref by paying a greater
manifold insertion ΔV (1 cm/s ×λ1 E10 m/s, as λ1 ¼1406 for Az ¼4×105 km). Hence, the minimum phasing ΔV under the relaxed time slot is the minimum among the three ΔVmin at the three θ0. Fig. 16 shows the ΔVmin as a function of θ0(i) . It is shown that the required ΔV to cover all θ0 becomes 35 m/s. In other words, for the most difficult θ0(i) (215° for Az ¼4×105 km), the spacecraft can depart in the next revolution with a ΔV saving of 12 m/s. In addition, it can be seen that, a half fraction of θ0(i) require a ΔV smaller than 5 m/s, while another half is a blackout region where the required ΔV is much greater. This characteristic can be considered to decide a launch date favorable for extended missions. The ΔV to cover all θ0 is theoretically zero if the manifold insertion timing is not constrained. To sum up, the minimum phasing ΔV vs θ0 as shown in Fig. 11 can be computed beforehand as a useful reference when planning the halo orbit mission and the extended missions that require a lunar encounter. Fig. 15. The optimal τ, t1 and tf.
6. Conclusions As several successful extended missions such as the ISEE-3 and
ARTEMIS applied transfers to the Moon, this paper investigates the minimum ΔV cost for going to the Moon from Sun–Earth halo orbits for various initial lunar phase θ0. Due to the solution
472
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Acknowledgements The research was carried out at ISAS/JAXA under the joint research program between JAXA and Japan's universities. The first author appreciates the resources provided by JAXA, especially the advice from Dr. Stefano Campagnola and Dr. Chit Hong Yam. The first author also appreciates the support from her current workplace, Technology and Engineering Center for Space Utilization, Chinese Academy of Sciences, for finalizing this paper.
Appendix A
Fig. 16. Minimum phasing ΔV vs the θ0 (upper); minimum phasing ΔV vs the middle θ0 (lower) given a time slot of two halo revolutions for phasing planning.
multiplicity of the TPBVP, the constrained optimizer that solves TPBVP implicitly can possibly not converge to the solutions of the lowest ΔV case. A TPBVP solver aided by a two-body Lambert solver is developed. Optimizing with the TPBVP solver that identifies short-way and long-way approaches recovers the two corresponding variations of minimum ΔV with θ0. Constrained optimization tools (e.g. fmincon and SNOPT) are convenient to use, but may be problematic for some specific problems. For similar problems, constrained optimization can be employed at first. If there are minima that are suspected to be missed, the specific optimization with the TPBVP solver can be used for cross-checking suspicious intervals. Although the scenario in this paper is transferring from Sun–Earth L2 halo orbits to the Moon, the developed TPBVP solver can be used to identify multiple trajectories connecting two-body and three-body realms in other interplanetary applications. In particular, there are certainly applications to go to the Moon from the Sun–Earth realm, as demonstrated by the missions HITEN and WIND. The required phasing ΔV to cover all initial lunar phases for a lunar encounter does not significantly change as Az changes from 3.5×105 km to 5×105 km, which is around 45 m/s. Down to the point Az ¼ 3.5×105 km, required ΔV decreases significantly as Az decreases. The required ΔV is only 14 m/s for Az ¼1×105 km. The order of the required phasing ΔV is considered acceptable for extended missions. Furthermore, the trajectories with minor phasing ΔV resemble the four types of ballistic manifold transfer trajectories to the Moon. These four types of transfer trajectories have different properties in lunar gravity assist and lunar mission scenarios. Therefore, one can choose from the four options depending on the condition and mission requirement. In particular, there are two transfer options applicable for lunar missions. It is computed that the sum of ΔVmin and V1 with respect to the Moon is around 500 m/s for the example halo orbit of an Az of 4×105 km. The paper also discusses the phasing planning based on the result of minimum ΔV and the shift of initial lunar phase with halo orbit revolution. The variation of the minimum phasing ΔV with lunar phase is a useful reference, when it comes to planning a Sun–Earth halo orbit mission and the extended mission involving a lunar encounter.
The eigenvalues of the monodromy matrix of a halo orbit consists of a pair of real eigenvalues λ1 and λ2 with λ1 λ2 ¼ 1, a pair of conjugate complex eigenvalues λ3 and λ4 with modulus 1, and a pair λ5 and λ6 equal to unity. The eigenvector e1 associated to the dominant eigenvalue λ1 (λ1 41) points to the divergent direction, and similarly e2 associated to λ2 points to the convergent direction. The eigenvectors e3 and e4 of the complex eigenvalues define a span formed by erotsum = (e3 + e4 ) /2 and erotdiff = (e3 − e4 ) /2i . A small displacement falls in this span will rotate arg (λ3) in the span every period of the halo orbit. Because the eigenvalue of 1 has only one eigenvector, i.e. e5 ¼ e6, the algebraic multiplicity of the monodromy matrix is two but the geometric multiplicity is one. The monodromy matrix is defective. To complete the eigenspace, the generalized eigenvectors are derived as follows, 2 el6 = Ker ⎡⎣ ( M − λ 6 I) ⎤⎦
(A.1)
Because the algebraic multiplicity is two, there are two solutions to Eq. (A.1). To have geometric meaning, the vector e˜ 6 that ^ and perpendicular to e5 is selected to lies in the span of the two e 6 ˜ be the sixth eigenvector. e6 points to nearby periodic orbits. The normalized eigenvectors e¯ i (i ¼ 1,…,6) for e1, e2, erotdiff , erotdiff , e5 and e˜ 6 form the eigenspace as follows,
E = ⎡⎣ e¯ 1 e¯ 2 e¯ 3 e¯ 4 e¯ 5 e¯ 6 ⎤⎦
(A.2)
An infinitesimal displacement of magnitude c1 along e¯ 1 will grow λ1 times in a period. Small displacements along other eigenvectors e¯ i (i¼ 2,…, 6) will not diverge as the modulus of the corresponding eigenvalues are smaller than or equal to 1. Then, as long as the components of the perturbation along each eigenvector is small,
X (X 0 + c1e¯ 1, t ) ≈ X (X 0 + EC, t )
(A.3)
where X0 is the initial state vector (manifold insertion point along the halo orbit), X is the state vector at time t, and C is the coefficient vector T C = ⎡⎣ c1 c2 c3 c4 c5 c6 ⎤⎦
(A.4)
Due to this feature, for the c1 of a given manifold, one can compute the coefficients c2,…,c6 to zero position perturbations and therefore find the Δv to approximate the effect due to sixdimensional perturbation c1e¯ 1. The constraint is expressed as,
⎡ 0(3 × 1) ⎤ ⎢ ⎥ = EC ⎣ Δv(3 × 1) ⎦
(A.5)
The matrix E is partitioned into four submatrices as,
⎡ E11 (3 × 3) E12 (3 × 3) ⎤ ⎥ E=⎢ ⎢⎣ E21 (3 × 3) E22 (3 × 3) ⎥⎦ Substituting Eq. (A.6) to Eq. (A.5) results in,
(A.6)
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
⎡ c4 ⎤ ⎡ c1 ⎤ 0 = E11⎢ c2 ⎥ + E12 ⎢ c5 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ c6 ⎦ ⎣ c3⎦
References
(A.7)
and
⎡ c4 ⎤ ⎡ c1 ⎤ Δv = E21⎢ c2 ⎥ + E22 ⎢ c5 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ c6 ⎦ ⎣ c3⎦
(A.8)
Integrating Eqs. (A.7) and (A.8) to cancel [c4 c5 c6]T yields the expression for Δv; that is,
⎡ c1 ⎤ Δv = A (3 × 3) ⎢ c2 ⎥ ⎢ ⎥ ⎣ c3⎦
(A.9)
where, −1 A = (E21 − E22 E12 E11)
(A.10)
c1 is associated to the desired the manifold, and the combination of c2 and c3 determines the magnitude of Δv. Let the entry in the ith row and j-th column of A be defined as aij. Based on Eq. (A.9), the square of ΔV can be expressed as,
Δv 2 = b1c12 + b2 c22 + b3 c32 + b12 c1c2 + b23 c2 c3 + b13 c1c3
(A.11)
where, 3
bj =
3
∑ aij bkj = 2 ∑ aik aij i=1
i=1
(A.12)
At the minimum, the derivatives of Eq. (A.11) with respect to c2 and c3 should be zero; that is,
∂ΔV 2/∂c2 = b12 c1 + 2b2 c2 + b23 c3 = 0
(A.13)
∂ΔV 2/∂c3 = b13 c1 + b23 c2 + 2b3 c3 = 0
(A.14)
Then, c2 and c3 can be calculated from, −1 ⎡ c2 ⎤ ⎡ 2b2 b23 ⎤ ⎡ − b12 ⎤ ⎥ ⎢ ⎥ c1 ⎢⎣ c ⎥⎦ = ⎢ 3 ⎣ b23 2b3⎦ ⎣ − b13 ⎦
(A.15)
Let B(2×1) represent the resulting vector before c1. Then, the minimum required Δv corresponding to the manifold of c1 is expressed as,
⎡1 ⎤ Δv = A ⎢ ⎥c1 ⎣ B (2 × 1) ⎦ Meanwhile, the optimal direction to place folds is,
⎡1 ⎤ Δl v = A⎢ ⎥ ⎣ B (2 × 1) ⎦
473
(A.16)
Δv to generate mani-
(A.17)
This ΔV direction corresponds to the 6-dimensional divergent direction of the halo orbit.
[1] R.W. Farquhar, The flight of ISEE-3/ICE : origins, mission history, and a legacy, in: Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 1998. doi: 〈http://dx.doi.org/10.2514/6.1998-4464〉. [2] E.A. Belbruno, J.K. Millert, Sun-perturbed Earth-to-Moon transfers with ballistic capture, J. Guid. Control. Dyn. 16 (1993) 770–775. [3] W.S. Koon, M.W. Lo, J.E. Marsden, S.D. Ross, Low energy transfer to the Moon, Dyn. Nat. Artif. Celest. Bodies (2001) 63–73, http://dx.doi.org/10.1007/ 978-94-017-1327-6_8. [4] K.C. Howell, M. Kakoi, Transfers between the Earth–Moon and Sun–Earth systems using manifolds and transit orbits, Acta Astronaut. 59 (2006) 367–380, http://dx.doi.org/10.1016/j.actaastro.2006.02.010. [5] G. Mingotti, F. Topputo, F. Bernelli-Zazzera, Optimal low-thrust invariant manifold trajectories via attainable sets, J. Guid. Control. Dyn. 34 (2011) 1644–1656, http://dx.doi.org/10.2514/1.52493. [6] C. Martin, B.A. Conway, Optimal low-thrust trajectories using stable manifolds, in: B.A. Conway (Ed.), Spacecraft Trajectory Optimization, 1st ed.,Cambridge University Press, 2010, pp. 238–262 http://dx.doi.org/10.1017/CBO9780511778025.010. [7] D.C. Folta, M. Woodard, K. Howell, C. Patterson, W. Schlei, Applications of multi-body dynamical environments: the ARTEMIS transfer trajectory design, Acta Astronaut. 73 (2012) 237–249, http://dx.doi.org/10.1016/j. actaastro.2011.11.007. [8] M. Nakamiya, H. Yamakawa, D.J. Scheeres, M. Yoshikawa, Interplanetary transfers between halo orbits: connectivity between escape and capture trajectories, J. Guid. Control. Dyn. 33 (2010) 803–813, http://dx.doi.org/10.2514/ 1.46446. [9] G. Mingotti, F. Topputo, F. Bernelli-Zazzera, Earth–Mars transfers with ballistic escape and low-thrust capture, Celest. Mech. Dyn. Astronaut. 110 (2011) 169–188, http://dx.doi.org/10.1007/s10569-011-9343-5. [10] J.J. Guzman, D.W. Dunham, P.J. Sharer, J.W. Hunt, C.J. Ray, H.S. Shapiro, et al., STEREO mission design implementation, in: Proceedings of the 20th International Symposium on Spacecraft Flight Dynamics, Annapolis, Maryland, USA, 2007. 〈http://issfd.org/ISSFD_2007/21-4.pdf〉. [11] J. Kawaguchi, H. Yamakawa, T. Uesugi, H. Matsuo, On making use of lunar and solar gravity assists in LUNAR-A, PLANET-B missions, Acta Astronaut. 35 (1995) 633–642, http://dx.doi.org/10.1016/0094-5765(95)00013-P. [12] S. Campagnola, C. Corral, R. Jehn, Design of lunar gravity assist for the bepicolombo mission to mercury, Adv. Astronaut. Sci. 119 (2004) 427–442. [13] D. Landau, T.P. Mcelrath, D. Grebow, N.J. Strange, Efficient lunar gravity assists for solar electric propulsion missions, Adv. Astronaut. Sci. 143 (2012) 917–934. [14] T.P. McElrath, G. Lantoine, D. Landau, D. Grebow, N. Strange, R. Wilson, et al., Using gravity assists in the Earth–Moon system as a gateway to the solar system, in: Proceedings of the Global Space Exploration Conference, Washington, D.C., 2012. [15] H. Chen, Y. Kawakatsu, T. Hanada, Earth Escape from a Sun-Earth halo orbit using unstable manifold and lunar swingbys, Trans. Jpn. Soc. Aeronaut. Space Sci., 59, 2016. [16] Y. Kawakatsu, Destiny mission overview: a small satellite mission for deep space exploration technology demonstration, Adv. Astronaut. Sci. 146 (2012) 727–739. [17] V. Szebehely, Theory of Orbits: The Restricted Problem of Three Bodies, Academic Press, New York, San Francisco, London 1967, pp. 134–139. [18] W.S. Koon, M.W. Lo, J.E. Marsden, S.D. Ross, Dynamical systems, Three-Body Probl. Sp. Mission Des. (2011) 174–175. [19] K.C. Howell, Three-dimensional, periodic, “Halo” orbits, Celest. Mech. (1984). [20] G. Lantoine, T.P. Mcelrath, Families of solar-perturbed Moon-to-Moon transfers, in: Proceedings of theAAS/AIAA Space Flight Mechanics Meeting, 2014. [21] R.R. Bate, D.D. Mueller, J.E. White, Fundamentals of astrodynamics, Dover (1971). [22] N. Arora, R.P. Russell, N.J. Strange, Partial Derivatives of the Lambert Problem, in: Proceedings of the AIAA/AAS Astrodyne Specialist Conference, American Institute of Aeronautics and Astronautics, 2014. 〈http://dx.doi.org/10.2514/6. 2014-4427〉. [23] H. Chen, NEW Results of Orbital Dynamics and the Application to Orbit Prediction and Mission Design, Kyushu University 2015, pp. 148–153. [24] F. Topputo, R. Zhang, F. Bernelli-Zazzera, J. Luo, Numerical approximation of invariant manifolds in the restricted three-body problem, in: Proceedings of the International Astronautical Congress, Beijing, China, 2013. [25] R. Keys, Cubic convolution interpolation for digital image processing, Acoust. Speech Signal Process. IEEE Trans 29 (1981) 1153–1160, http://dx.doi.org/ 10.1109/TASSP.1981.1163711.
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/aa
Phasing Delta-V for transfers from Sun–Earth halo orbits to the Moon Hongru Chen a,n, Yasuhiro Kawakatsu b, Toshiya Hanada a a b
Department of Aeronautics and Astronautics, Kyushu University, Fukuoka, Japan Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency (ISAS/JAXA), Sagamihara, Japan
art ic l e i nf o
a b s t r a c t
Article history: Received 8 June 2015 Received in revised form 9 December 2015 Accepted 3 May 2016 Available online 26 May 2016
Inspired by successful extended missions such as the ISEE-3, an investigation for the extended mission that involves a lunar encounter following a Sun-Earth halo orbit mission is considered valuable. Most previous studies present the orbit-to-orbit transfers where the lunar phase is not considered. Intended for extended missions, the present work aims to solve for the minimum phasing ΔV for various initial lunar phases. Due to the solution multiplicity of the two-point boundary value problem, the general constrained optimization algorithm that does not identify multiple feasible solutions is shown to miss minima. A two-step differential corrector with a two-body Lambert solver is developed for identifying multiple solutions. The minimum ΔV associated with the short-way and long-way approaches can be recovered. It is acquired that the required ΔV to cover all initial lunar phases is around 45 m/s for the halo orbit with out-of-plane amplitude Az greater than 3.5×105 km, and 14 m/s for a small halo orbit with Az ¼ 1×105 km. In addition, the paper discusses the phasing planning based on the ΔV result and the shift of lunar phase with halo orbit revolution. & 2016 IAA. Published by Elsevier Ltd. All rights reserved.
Keywords: Extended missions Circular restricted three-body problem Halo orbits Lunar transfers Two-point boundary value problem
1. Introduction The halo orbit about a libration point in a three-body system has features such as relatively constant distances and orientation with respect to the two primary bodies, which are favorable for scientific observation, thermal control and communication. A famous mission is the ISEE-3. ISEE-3 achieved several missions including monitoring solar wind in a Sun–Earth halo orbit, exploring geo-tail and visiting comets beyond original expectations [1]. The success of the mission is owing to the elaborate planning of transfers where the three-body dynamics and multiple lunar gravity assists were used. There have been several libration point missions in addition to the ISEE-3, and many are being planned. The success of ISEE-3 and the increasing interest in halo orbit missions give rise to the need of investigation of extended missions following them. Utilizing the low-energy transfer along the manifold associated with a libration point orbit, various missions can be developed. One option for the extended mission of a Sun–Earth halo orbit mission is going to the Moon. Following the success of the HITEN mission, several studies and mission design methods have demonstrated the low-energy transfers to the Moon under the gravitational influences of the Sun and the Moon [2–6]. The extended mission of THEMIS, ARTEMIS, reached to the Sun–Earth L1 n
Corresponding author. E-mail address: [email protected] (H. Chen).
http://dx.doi.org/10.1016/j.actaastro.2016.05.003 0094-5765/& 2016 IAA. Published by Elsevier Ltd. All rights reserved.
region and then an Earth-Moon L2 Lissajous orbit applying manifold dynamics [7]. Another promising option is the interplanetary exploration. However, low-energy escape only permits accessing a few asteroids in a narrow zone. For instance, the low-energy transfer from the Earth to the Mars needs an additional ΔV of around 2 km/s [8,9]. On the other hand, lunar gravity assists have proven to be an effective way to increase escape energy [10–14]. The ISEE-3 transferred from the Sun–Earth L1 halo orbit to the Moon and then performed multiple lunar gravity assists for Earth escape. A previous study investigates the combination of the unstable manifold and multiple lunar swingbys for Earth escape [15]. The key remarks obtained from that study are: 1) manifold-guided lunar swingbys can achieve much higher characteristic energy (C3) with respect to the Earth than direct escape along the manifold; 2) for the case that the V1 with respect to the Moon is not sufficient for high-energy Earth escape, a second lunar swingby can efficiently increase the C3 to a level of 2.7 km2/s2 in 70–90 days. Transfers to the Moon utilizing manifold dynamics have been intensively studied and well demonstrated [2–6]. However, previous studies are to find low-cost transfer paths and meanwhile decide the optimal departure orbit phase or terminal lunar phase. Using numerical techniques, Ref. [4] solves the zero-ΔV transfer from a Sun–Earth L2 halo orbit to an Earth–Moon L1 halo orbit. Different from those studies on the orbit-to-orbit transfers, the motivation of this paper is the extended mission where the departure orbit, i.e. a Sun–Earth halo orbit, is not pre-phased for a lunar encounter. The previous work (Ref. [15]) assumes the transfers to the Moon are also ballistic. The present paper is a
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
follow-up work concerning a phasing maneuver for encountering the Moon. Some phasing ΔV were used to target the Moon in ISEE3 and ARTEMIS. Obtaining the required phasing ΔV for different conditions will facilitate the design of similar extended missions that involve a lunar encounter. The study is also motivated by JAXA's mission DESTINY (Demonstration and Experiment of Space Technology for INterplanetary voYage) [16]. DESTINY will be launched to a halo orbit about the Sun–Earth L2 point in 2021, and the possibility of extending the mission to visit a heliocentric body is currently under discussion. Section 2 presents the basis of the circular restricted threebody problem. Section 3 gives the overview and parameterization of the phasing problem for the transfer from a Sun–Earth halo orbit to the Moon. Section 4 presents the computation procedures and the results of minimum required phasing ΔV for various initial lunar phases. In attempts of optimization, the problem of missing minima is found due to the solution multiplicity of the two-point boundary value problem (TPBVP). To tackle this problem, a sequential differential corrector aided by a two-body Lambert solver is developed to identify the multiple solutions to a TPBVP robustly. The required ΔV to cover all initial lunar phases is acquired. In addition, the characteristics of the transfer trajectories with minor phasing ΔV are discussed for different mission scenarios, namely, lunar mission and lunar gravity-assisted Earth escape. Section 5 discusses the phasing planning based on the result of minimum ΔV, and the relationship between the initial lunar phase and the halo orbit revolution. Conclusions are given in the final section.
2. System dynamics 2.1. Circular restricted three-body problem Trajectories are considered in the circular restricted three-body problem (CR3BP) of the Sun–Earth system. The CR3BP assumes two primary bodies m1 and m2 moving in a circular orbit about their barycenter. The mass of the third body is negligible compared to the masses of the two primaries. The two primaries can be the Sun and Earth, the Earth and Moon or etc. The rotating coordinate system with the barycenter at the origin and the primaries fixed on the x-axis is chosen to describe the motion of the third body. For convenience, the angular velocity of the rotating frame, the total mass and the distance between the two primaries are normalized to 1. μ is the ratio of the mass of m2 to the mass of m1. Then, masses of m1 and m2 become 1-μ and μ in the normalized model. The coordinates of m1 and m2 are [ μ, 0, 0] and [1 μ, 0, 0], respectively, as shown in Fig. 1. The equations of motion of the third body are:
x¨ − 2ẏ = ∂U /∂x
(1)
y¨ + 2ẋ = ∂U /∂y
(2)
z¨ = ∂U /∂z
(3)
where the gravitational potential U in the system is
U = (x2 + y2 )/2 + (1 − μ)/rP1 + μ /r P2
(4)
where
r P1 =
(x + μ)2 + y2 + z2 , r P2 =
(x − 1 + μ)2 + y2 + z2
(5)
465
Fig. 1. Circular restricted three-body system and Lagrangian points.
from the equations of motion. They are also referred to as the libration points or Lagrangian points labeled as L1,…, L5. The three collinear points L1, L2 and L3 are unstable. The geometries of the five equilibrium points are depicted in Fig. 1. 2.2. Halo orbits and invariant manifolds The linearized equations at the collinear equilibrium points reveal that the motion in the vicinity of these points consists of the in-plane and out-of-plane oscillatory modes. A halo orbit is a periodic orbit with a sufficiently large out-of-plane amplitude Az such that the out-of-plane oscillatory frequency becomes equal to the in-plane oscillatory frequency [18]. Halo orbits are generally computed numerically [19]. The monodromy matrix is the state transition matrix for one period of a periodic orbit. By examining the eigenvalues and eigenvectors of the monodromy matrix, the insight of the orbit dynamics can be revealed. For halo orbits, there is a pair of eigenvalues with λ1 λ2 ¼ 1 [18]. The dominant eigenvalue λ1 is around 1500 for Sun–Earth L1/L2 halo orbits. A small displacement along the eigenvector of the dominant eigenvalue (divergent eigenvector for short) will propagate into exponentially increasing divergence from the halo orbit quickly, since the displacement will grow λ1 times in a period. Such a departure trajectory follows the unstable invariant manifold. A derived manifold trajectory can go earthward or anti-earthward. On the other hand, impulsive changes of both position and velocity are not practical in real missions. Δv are added to generate manifold trajectories instead. The Δv direction corresponding to the divergent eigenvector is computed based on Appendix A. To globalize the manifold, small Δv in the divergent directions are added to a series of points along the halo orbit. The propagation of the perturbed states yields an unstable manifold with a tube structure, as is shown in Fig. 2. Conversely, backward propagating the states perturbed by small Δv in convergent directions corresponding to the eigenvalue smaller than 1 (λ2) will result in the stable manifold that asymptotically converges into the halo orbit.
3. Phasing problem
An energy quantity,
J = 2U − (x2̇ + y2̇ + z2̇ )
(6)
holds in the system. It is known as the Jacobi integral. Ref. [17] presents the detail of resolving the equilibrium points
As both of the unstable invariant manifold and lunar gravity assist are useful in delivering a spacecraft to deep space, our previous study applies the combination of them for Earth escape. That study shows this strategy can achieve a considerable C3 for
466
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Az = 4x105 km
−3
6
x 10
corrected trajectory ωΔtf
Δxf, Δtf
rM (tf) xf, tf
4 1
initial trajectory (along manifold)
y ,AU
2 0
2
L2
4
Moon’s orbit
3
x 1 -, t 1 -
x1+, t1+ Δv
−2 Moon’s orbit
−4 −6
Fig. 3. Correct a manifold trajectory to target the Moon.
0.995
1
1.005
1.01
x ,AU −3
x 10
z, AU
2 2
0
3
1
L2
4
−2 0.995
1
1.005
1.01
x, AU Fig. 2. The intersection of the manifold tube associated with the halo orbit with Az ¼ 4×105 km and Moon's orbit exhibits four types of intersection.
interplanetary missions. An overview of the previous study is presented in this section, which leads to the discussion on the phasing problem. 3.1. Background The spacecraft is assumed initially in a Sun–Earth L2 halo orbit. The inclination and eccentricity of the Moon's orbit are small and assumed zero. This simplified model is sufficient for acquiring reliable information for practical mission design [20]. Fig. 2 shows the earthward trajectories of the unstable manifold of the halo orbits with Az ¼4×105 km and the Moon's orbit. The grey circles mark the intersections of the manifold trajectories and the ecliptic plane, outlining the intersection lines of the manifold tube and the ecliptic plane. It is found that the manifold and the Moon's orbit intersect at four points after the first crossing through Moon’s orbit. These four types of intersection exist for the wide Az range investigated (from 1×105 km to 5×105 km). The four types of intersection are numbered in order of the lunar phases. The previous study shows, if the Moon is exactly at the intersection points at the manifold arrival, a following single or double lunar swingby can lift the C3 energy with respect to the Earth to around 2.7 km2/s2, which is far greater than a direct escape from the halo orbit via the anti-earthward unstable manifold. However, since the manifold tube is fixed for a specified halo orbit, the positions of the four intersection points are fixed. The Moon may not be at the intersection points at the time the manifold reaches. Therefore, a phasing maneuver to achieve a lunar encounter is considered, and the cost of it is investigated in this paper. 3.2. Parameterization Define the epoch the spacecraft passes a reference point rref in the halo orbit as the reference epoch tref. rref is set at [x L2 + Ax , 0, Az ] in this study, where [x L2, 0, 0] is the coordinate of L2
and Ax is the x-amplitude of the orbit. Define the lunar phase at tref as the initial lunar phase θ0. θ0 is related with the launch date and halo orbit revolution. The halo orbit size Az, which is associated with the Jacobi integral J is known once the mission is decided. For a given Az, the required ΔV to target the Moon for various θ0 would be a useful reference for evaluating and planning the extended mission involving a lunar encounter. It is reasonable to assume that the scheme of leaving the halo orbit via the manifold and then applying a mid-course ΔV to target the Moon is an optimal choice. To find the optimal patch point along the manifold, the manifold is parameterized. Two time variables τ and t1 are introduced. τ is the time of the manifold insertion point from tref. The magnitude of manifold insertion Δv is set to be 1 cm/s, then, τ determines a manifold trajectory. t1 is the propagation time of the manifold and also the timing of the phasing maneuver. The final epoch tf determines the lunar encounter position rM, which is a function of θ0 and τ þtf. The maneuver at t1 corrects the manifold trajectory to reach rM at tf. Fig. 3 schematically depicts the phasing problem. A control of V1 with respect to the Moon can further enhance the mission design. The V1 is a critical factor to the mission using lunar gravity assists. However, for a lunar mission, the V1 should be reduced. The control of V1 should depend on the specific mission scenario. Discussions for various mission scenarios are beyond the scope of this paper. As for using lunar gravity assists to achieve Earth escape, studies have shown that even if the V1 is small, a second lunar swingby without additional ΔV cost can efficiently lift the escape energy to the theoretical maximum level in 2-3 months [14]. Therefore, the objective of phasing in this study is simply to achieve a lunar encounter with minimum ΔV. However, characteristics of the lunar transfer trajectories with minor phasing ΔV will be revealed and discussed for different mission scenarios in the subsequent section.
4. Optimization of phasing ΔV 4.1. Constrained optimization Trajectories are subject to Eqs. (1)–(3). For given Az and optimization is to find the 6-element set, T y = ⎡⎣ τ , t1, tf , ΔvT⎤⎦
θ0, the (7)
that satisfies the equality boundary constraint,
c = r M (tf ) − rf = 0
(8)
where rf is the position of the spacecraft at tf, and at the same time minimizes the cost function,
ΔV = Δv
(9)
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
467
A number of initial guesses of y should be tried in order to find the global minimum. It is prescribed that,
−T < τ < 0
(10)
where T is the period of the halo orbit, which is a function of Az. This prescription suggests a slot of one halo orbit period before tref for targeting a lunar encounter. 20 τ′ are sampled evenly within the prescribed slot as the initial guesses. (Prime and asterisk superscripts are used to indicate initial guesses and optimal solutions, respectively.) Then, the full range of manifold trajectories is taken into account. Note that, once τ* is solved, it does not mean the spacecraft must be inserted to the manifold at τ*. If the spacecraft misses τ*, it can be inserted to the same manifold trajectory at a later τ by paying a correspondingly larger manifold insertion ΔV. As mentioned in Section 3.1, encounter chances and closest approaches can occur at every intersection of the manifold tube and the Moon's orbit. There can be infinite local minima if tf is not constrained. To simplify the problem, this study only concerns the first encounter chances. Then, tf′ is chosen to be the epoch of the first closest approach to the Moon along the manifold trajectory. As for t1′, through experiments with some boundary conditions, it is observed that there are usually two local minima of ΔV along the t1 coordinate (an example will be presented in Fig. 7). t1′ are chosen to be 30 and 120 days prior to tf′ based on the observation. It becomes difficult for the numerical computation to converge as t1′ approaches 0, because of the high sensitivity in the divergent direction in the halo orbit regime (λ1 E 1500 compared to other eigenvalues r 1). Therefore, the ranges and sampled initial guesses of τ, t1 and tf are considered reasonable and sufficient for the problem setting. Δv′ is 0. The MATLAB constrained optimization solver, fmincon, is used in this study. Minimum ΔV is computed for the θ0 varied from 0° to 360°. θ0 is sampled every 5°. Eventually, the required ΔV budget to cover the full range of θ0 can be obtained, providing a reference for extended missions. Fig. 4 summarizes the procedures to generate the minimum phasing ΔV as a function of θ0. 4.2. Results and discussions of constrained optimization The halo orbit with Az ¼4×105 km is taken as an example to investigate the minimum phasing ΔV, ΔVmin, as a function of θ0.
Fig. 5. The minimum phasing ΔV and the optimal τ solved by the constrained optimization (Az ¼ 4×105 km).
The ΔVmin vs θ0 computed by the workflow in Fig. 4 is shown in Fig. 5. The ΔVmin profile is shown made up of several profiles. Let τ1,…, τ4 represent the τ of the four manifold trajectories that intersect with the Moon's orbit, and θ01…, θ04 correspondingly represent the θ0 conditions under which the four manifold trajectories can meet the Moon without phasing ΔV. The τ* of each profile are found near τ1,…, τ4 for the θ0 near θ01…, θ04 (see the τ* variation measured by the scale to the right), respectively. Then, the minima can be grouped into four groups corresponding to the four types of intersection. To get an insight, Fig. 6 shows the ΔVmin variations associated with Type 1 and Type 2 in the θ0 interval [300° 340°]. It can be seen that the ΔVmin associated with the two types take turns to stand out (i.e., to have lower ΔV for the θ0 by comparison) in this interval. Because a sufficient number of τ′ are tried in the computation and the τ of the two types are quite different, the program can successfully identify the minima of the two groups and output the lower ones. On the other hand, it is observed that the zero minimum of the Type-3 intersection as well as the minima associated with it are missing in the result. This problem is found related with the solution multiplicity of the two-point boundary value problem (TPBVP). It is well known that, in the two-body problem, there are multiple solutions to the TPBVP, namely, the short-way, long-way and multiple-revolution approaches. The manifold trajectory of the Type-3 intersection can be regarded as a long-way approach in the near-Earth realm; while the other three are short-way approaches (see Fig. 2). For the group-3 and group-4 manifold trajectories, the required phasing ΔV that leads to short-way and long-way approaches may not differ greatly. For illustration, Fig. 7 shows the ΔV variations associated with the two approaches for various two points, which are defined by sets of t1 and tf, for a given θ0 near the θ03 and a given τ near τ3. It can be seen that, each solution profile exhibits two local minima. Moreover, for a considerable range, the ΔV for the two approaches are close to each other. The constrained optimization algorithm solves TPBVP implicitly. Once the optimizer converges to a feasible solution, it is inevitable for the gradient-based approach to be trapped in the current solution regime. Therefore, the lowest ΔV minimum can be possibly missed by a general constrained optimizer. One may try multiple initial guesses for Δv. However, as there are actually three variables in Δv, if no good knowledge on their ranges is available, it may take excessively long time on many random initial guesses in order to guarantee all feasible solutions. Good estimates of the Δv of both trajectory motions are desirable in this problem. A TPBVP solver is developed in the next subsection for the purpose of distinguishing different solutions. 4.3. Sequential differential corrector with a Lambert solver for TPBVP
Fig. 4. Workflow of computing the minimum phasing ΔV vs θ0 for a given Az using constrained optimization.
To identify the Δv for the desired trajectory motion, the first step is to obtain a close guess utilizing a two-body Lambert solver
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Group 1 Group 2 (near τ2 = -144 days) (near τ1 = -17.7 days)
468
ΔVmin(group 1) > ΔVmin(group 2)
ΔVmin(group 1) < ΔVmin(group 2)
ΔVmin(group 1) > ΔVmin(group 2)
Fig. 6. The minimum ΔV of group 1 and group 2 for the θ0 varied from 300° to 340° (Blank areas are where the ΔV is too high to be displayed or failed to be solved by the program).
Fig. 7. ΔV contours of the short-way and long-way approaches. (Blank areas are where ΔV is too high to be displayed or failed to be solved by the program).
that can identify multiple solutions in the two-body problem. The trajectory is segmented at an auxiliary epoch t2 when the manifold trajectory enters a defined Earth-centered sphere with a radius of 8×105 km. The sphere is considered as the boundary between the three-body and two-body realms. Then, the t1 t2 arc is seen as a three-body trajectory, while the t2 tf arc a two-body trajectory. The t1 t2 arc is propagated, which gives the position r2 and velocity v2 at t2. Let "i" superscript indicate vectors referenced to the inertial frame. The t2 tf arc connecting r2 and rfi is solved by a two-body Lambert solver. rfi is the rf rotated by an angle of Ω (tf t2), where Ω is the rotating rate of the Sun Earth system. The Lambert solver used in this study is based on Bate's method [21] and can identify the solution (v2i+) for a desired (short-way or long-way) trajectory motion. Then, the corresponding v2+ referenced to rotating frame can be obtained. There appears an inconsistency,
Fig. 8. Break the trajectory at an auxiliary epoch t2, and then smoothen the connection of the two-body and three-body arcs by differential correction.
Δv2 = v 2+ − v 2−
(11)
By reducing Δv2, a good initial guess for the Δv to establish the desired trajectory can be obtained. The v2 can be tuned by a differential correction based on the propagated state transition matrix of the t1 t2 arc. The method of deriving the state transition matrix for CR3BP can be found in Ref. [19]. A recent study [22] gives the analytic expressions of the partial derivatives of the twobody Lambert problem, providing the gradients to correct v2 þ as well. Fig. 8 illustrates the correction concept. In detail, for the t1 t2 arc, differentiating the state x2 at t2 with respect to the state x1 at t1 yields,
δ x2 = Φ (t2, t1) δ x1
(12)
where Φ(t2, t1) is the 6×6 state transition matrix that maps the state at t1 to the state at t2. Eq. (12) can be expanded to
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
⎡ δr2 ⎤ ⎡ Φ rr (3 × 3) Φ rv (3 × 3) ⎤ ⎡ δr1 ⎤ ⎥⎢ ⎥ ⎢ ⎥=⎢ ⎣ δ v2 ⎦ ⎣ Φ vr (3 × 3) Φ vv (3 × 3) ⎦ ⎣ δ v1⎦
−3
x 10
(13)
5
((t2, t1) is dropped for brevity). Then, the partial derivative of v2 with respect to Δv is
4
∂v −2 = Φ vv ∂Δv
2
(15)
The r2 change results in two folds of v2 þ changes. One is the change of inertial velocity V 2i + computed by the Lambert solver. The other is the change of the frame-induced item Δv i2→ r ( r2), which is added to the inertial velocity to acquire the velocity referenced to the rotating frame; that is,
v 2+ = v i2+ + Δv 2i→ r
(16)
Δv i2→ r ( r2) is calculated from
Δv 2i→ r
⎡ 0 1 0⎤ ⎢ ⎥ = ⎢ − 1 0 0⎥ r2 = Rr2 ⎣0 0 0⎦
(17)
(18)
The complete derivation and expression for L are given in Ref. [22]. With Eqs. (15), (17) and (18), the derivative of Eq. (16) with respect to Δv can be written as,
∂v 2+ = (L + R) Φ rv ∂Δv
(19)
Then, integrating Eqs. (11), (14) and (19) results in the correction for Δv at the i-th iteration, which is
Δv (i) − Δv (i − 1) = (Φ vv − (L + R) Φ rv )−1Δv 2(i − 1)
y, AU
1 tf Moon
0 −1
short way
long way
t1
−2 −3 −4 −5 0.998
1
1.002
1.004
1.006
1.008
1.01
x, AU
Fig. 9. An example of the short-way and long-way approaches to the Moon identified by the TPBVP solver.
4.4. Unconstrained optimization with the developed TPBVP solver
A partial derivative of the Lambert problem is,
L = ∂v i2+/∂r2
3
(14)
The change of Δv will result in a change of r2, which is expressed as,
∂r2 = Φ rv ∂Δv
469
With the TPBVP solver developed to identify different trajectory motions, the boundary condition is always satisfied. Hence, the unconstrained optimization can be applied. The MATLAB unconstrained optimizer, fminunc, is used. As Δv is solved by the TPBVP solver, the number of the independent variables becomes three, which are, T y = ⎡⎣ τ , t1, tf ⎤⎦
(22)
Fig. 10 shows the procedures of using the developed TPBVP solver and unconstrained optimization. An effort was made to derive the first-order derivatives (G¼ ∂Δv/∂y) and second-order derivatives (H ¼∂2Δv/∂y2) to speed up the computation and ensure accuracy. The expressions of the first- and second-order derivatives are
(20)
Once Δv2 meets the tolerance (30 m/s in this study), Δv is considered a close estimate of the true solution. Note that there are situations that Δv2 cannot be reduced to the tolerance in many (100 is used in this work) iterations. A failure to converge suggests a large discontinuity between v2 and v2 þ , as well as a large ΔV for the desired trajectory. As the main interest is the minimum ΔV for given Az and θ0, for this case, the first step will assign an infinite value to the ΔV for the current guess, without going to the next step of correction. In addition, in the process of optimization, t1 may be tuned to be smaller than t2. In that case, a good Δv estimate is directly obtained from the two-body Lambert solver. The second step is correcting the Δv estimate obtained from the previous step, to acquire the accurate Δv for the three-body t1 tf trajectory. The trajectory with the estimated Δv is propagated to tf, and differential correction is carried out iteratively to gradually zero the position difference at tf (Eq. (8)) based on Φ(tf, t1); that is,
Δv (i) − Δv (i − 1) = − Φ rv (t f , t1)−1c
(21)
With the Δv guess obtained from the first step aided by the Lambert solver, the differential correction on the three-body trajectory will robustly converge to the Δv for the desired trajectory motion. Fig. 9 shows an example of the short-way and long-way approaches identified for a set of two points.
Fig. 10. Workflow of computing the minimum phasing ΔV vs θ0 for a given Az using unconstrained optimization with the TPBVP solver identifying short-way or longway solution.
470
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Fig. 11. The minimum phasing ΔV vs θ0 and the optimal τ solved by the three optimization programs (Az ¼4×105 km).
presented in Ref. [23]. On the other hand, it is time-consuming to tune the Δv starting from 0 for every y determined by the optimizer. After the Δv for a given set of two points, which are defined by y(i 1), is solved, the initial guess of the Δv for the next neighboring y(i) (if ||y(i) y(i 1)|| o1 day) is given by,
Δv (i) ′ = G (i − 1) (y (i) − y (i − 1) ) + Δv (i − 1)
(23)
If y(i) is not regarded as a neighboring point to the previous one, Δv′ starts from 0. In addition, to avoid repeated time-consuming integration to τ and t1, the manifold states are approximated by fitting the surface to the states on the grid of sampled τ and t1 using the bicubic convolution interpolation [24,25].The convolution interpolation function can also give derivatives with respect to τ and t1, with which G and H can be computed [23]. 4.5. Results and discussions of unconstrained optimization Fig. 11 displays the ΔVmin solved by the constrained optimization, unconstrained optimization with short-way and long-way solution controls, and the corresponding τ*. Optimizing with the trajectory motion control recovers the two corresponding ΔVmin variations. The results of constrained optimization generally follow the short-way solutions. The optimal short-way solutions are shown associated with the intersections of Type 1, Type 2 and Type 4, while the optimal long-way solutions are associated with the intersection of Type 3. Integrating the two ΔV profiles, it can be seen that a ΔV budget of 47 m/s can cover all initial lunar phases for a lunar encounter for Az ¼ 4×105 km. Fig. 12 shows the required ΔV to cover all θ0 as a
function of Az. The required ΔV decreases as Az decreases. In particular, the required ΔV decreases relatively quickly with Az starting from Az ¼ 3.5×105 km. By inspection, the profile of group 1 moves leftward, and that of group 2 moves rightward as Az decreases, while the profiles of group 3 and group 4 do not vary much. Group 1 starts to stand out for the most difficult θ0 instead of group 2 at Az ¼3.5×105 km. As ΔVmin peaks occur at the intersections of the ΔVmin profiles, the leftward movement of the profile of group 1 leads to a faster decreasing rate of the peak ΔVmin. For a small halo orbit with Az ¼1×105 km, a ΔV of 14 m/s is sufficient for targeting the Moon. The order of the ΔV cost is considered acceptable. For instance, ISEE-3 processed a ΔV capacity of 310 m/s after completing the original mission, and a ΔV of 72 m/s was allocated to carry out the transfer from the halo orbit to the escape trajectory [1]. Ref. [15] has presented the properties of the four types of ballistic lunar encounter for Earth escape. The V1 of the intersections of Type 3 and Type 4 with respect to the Moon is around 1.35 km2/ s2, which is much higher than that of Type 1 and Type 2. The corresponding lunar encounters of Type 3 and Type 4 can lead to effective lunar gravity-assisted Earth escape. Although the V1 of Type 1 and Type 2 are only around 0.4 km2/s2, a second lunar swingby can efficiently lift up the escape energy with some extra flight time. In addition, the manifold trajectories of Type 1 and Type 2 resemble the low-energy transfers to the Moon presented in Refs. [2–4]. Therefore, for an extended lunar mission there are two applicable options. One can solve for the ΔVmin vs θ0 associated with the Type-1 and Type-2 intersections, by optimizing near τ1 and τ2. Fig. 13 shows the results. Fig. 14 shows the corresponding phasing trajectories for various θ0. As the period of the halo orbit is roughly six times the lunar period, for a certain type of intersection, a corresponding fraction of manifold trajectories can be targeted to the Moon with minor ΔV (o100 m/s). A minor ΔV does not change the shape of a trajectory and the final lunar encounter situation too much, but mainly alter the time of flight for waiting for or catching up with the Moon. Fig. 13 also displays the corresponding V1 with respect to the Moon and ΔVmin þV1, which can serve as an estimate of the cost of a lunar mission. If only phasing ΔV is optimized, with the two options of transferring to the Moon, ΔVmin þ V1 varies from 475 m/s to 550 m/s for different θ0. The Moon's gravity might contribute to reducing the required ΔV in a more realistic model. However, computing the ΔVmin for Type-1 and Type-2 intersections can give a quick estimate of the required ΔV in preliminary lunar mission design. Fig. 15 displays the variations of τ*, t1* and tf*. τ*, t1* and tf* evolve
Fig. 12. The required ΔV to cover all θ0 as a function of Az.
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
x 10
471
Type 1
-3
6 manifold departure
4 phasing trajectories
Moon y, AU
2
L2
0 2 Δv
4 6 1
x 10
1.01
Type 2
-3
6
Fig. 13. The ΔVmin vs θ0 associated with the Type-1 and Type-2 lunar intersections, the resulting V1 with respect to the Moon and ΔVmin þ V1 (Az ¼ 4×105 km).
phasing trajectories
4
Δv
2 y, AU
linearly with θ0. The above characteristics result in the linear change of ΔVmin with θ0, as is shown in Figs. 11 and 13.
1.005 x, AU
L2
0 Moon
5. Phasing planning 2
If the manifold insertion slot in practical missions is planned in previous or next halo revolutions, the obtained ΔV vs θ0 is still (i ) refer to the epoch of rref of the i-th halo orbit reuseful. Let tref volution. The period of the halo orbit (Az ¼4×105 km) is 179.7 days and the period of the Moon's orbit (with respect to the Sun–Earth rotating frame) is 29.5 days. It can be inferred that the θ0 shifts 33° (i ) at every tref . Therefore, one can relate θ0 to different tref , and find the halo orbit revolution and corresponding θ0 that requires acceptably low phasing ΔV. For instance, supposing that a two-revolution slot is permitted for inserting the spacecraft to the (i − 1) (i ) , tref and manifold, the lunar phases at the consecutive epochs tref
4
manifold departure
6 1
1.005
1.01
x, AU Fig. 14. The transfer trajectories with minimum phasing ΔV for the Type-1 and Type-2 intersections. Arrows represent phasing Δv.
(i + 1) are θ0(i − 1) , θ0(i) and θ0(i + 1) , respectively. It should be noted that tref
the optimal manifold trajectory solved for θ0(i − 1) can be targeted (i − 1) (i ) during the revolution between tref and tref by paying a greater
manifold insertion ΔV (1 cm/s ×λ1 E10 m/s, as λ1 ¼1406 for Az ¼4×105 km). Hence, the minimum phasing ΔV under the relaxed time slot is the minimum among the three ΔVmin at the three θ0. Fig. 16 shows the ΔVmin as a function of θ0(i) . It is shown that the required ΔV to cover all θ0 becomes 35 m/s. In other words, for the most difficult θ0(i) (215° for Az ¼4×105 km), the spacecraft can depart in the next revolution with a ΔV saving of 12 m/s. In addition, it can be seen that, a half fraction of θ0(i) require a ΔV smaller than 5 m/s, while another half is a blackout region where the required ΔV is much greater. This characteristic can be considered to decide a launch date favorable for extended missions. The ΔV to cover all θ0 is theoretically zero if the manifold insertion timing is not constrained. To sum up, the minimum phasing ΔV vs θ0 as shown in Fig. 11 can be computed beforehand as a useful reference when planning the halo orbit mission and the extended missions that require a lunar encounter. Fig. 15. The optimal τ, t1 and tf.
6. Conclusions As several successful extended missions such as the ISEE-3 and
ARTEMIS applied transfers to the Moon, this paper investigates the minimum ΔV cost for going to the Moon from Sun–Earth halo orbits for various initial lunar phase θ0. Due to the solution
472
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
Acknowledgements The research was carried out at ISAS/JAXA under the joint research program between JAXA and Japan's universities. The first author appreciates the resources provided by JAXA, especially the advice from Dr. Stefano Campagnola and Dr. Chit Hong Yam. The first author also appreciates the support from her current workplace, Technology and Engineering Center for Space Utilization, Chinese Academy of Sciences, for finalizing this paper.
Appendix A
Fig. 16. Minimum phasing ΔV vs the θ0 (upper); minimum phasing ΔV vs the middle θ0 (lower) given a time slot of two halo revolutions for phasing planning.
multiplicity of the TPBVP, the constrained optimizer that solves TPBVP implicitly can possibly not converge to the solutions of the lowest ΔV case. A TPBVP solver aided by a two-body Lambert solver is developed. Optimizing with the TPBVP solver that identifies short-way and long-way approaches recovers the two corresponding variations of minimum ΔV with θ0. Constrained optimization tools (e.g. fmincon and SNOPT) are convenient to use, but may be problematic for some specific problems. For similar problems, constrained optimization can be employed at first. If there are minima that are suspected to be missed, the specific optimization with the TPBVP solver can be used for cross-checking suspicious intervals. Although the scenario in this paper is transferring from Sun–Earth L2 halo orbits to the Moon, the developed TPBVP solver can be used to identify multiple trajectories connecting two-body and three-body realms in other interplanetary applications. In particular, there are certainly applications to go to the Moon from the Sun–Earth realm, as demonstrated by the missions HITEN and WIND. The required phasing ΔV to cover all initial lunar phases for a lunar encounter does not significantly change as Az changes from 3.5×105 km to 5×105 km, which is around 45 m/s. Down to the point Az ¼ 3.5×105 km, required ΔV decreases significantly as Az decreases. The required ΔV is only 14 m/s for Az ¼1×105 km. The order of the required phasing ΔV is considered acceptable for extended missions. Furthermore, the trajectories with minor phasing ΔV resemble the four types of ballistic manifold transfer trajectories to the Moon. These four types of transfer trajectories have different properties in lunar gravity assist and lunar mission scenarios. Therefore, one can choose from the four options depending on the condition and mission requirement. In particular, there are two transfer options applicable for lunar missions. It is computed that the sum of ΔVmin and V1 with respect to the Moon is around 500 m/s for the example halo orbit of an Az of 4×105 km. The paper also discusses the phasing planning based on the result of minimum ΔV and the shift of initial lunar phase with halo orbit revolution. The variation of the minimum phasing ΔV with lunar phase is a useful reference, when it comes to planning a Sun–Earth halo orbit mission and the extended mission involving a lunar encounter.
The eigenvalues of the monodromy matrix of a halo orbit consists of a pair of real eigenvalues λ1 and λ2 with λ1 λ2 ¼ 1, a pair of conjugate complex eigenvalues λ3 and λ4 with modulus 1, and a pair λ5 and λ6 equal to unity. The eigenvector e1 associated to the dominant eigenvalue λ1 (λ1 41) points to the divergent direction, and similarly e2 associated to λ2 points to the convergent direction. The eigenvectors e3 and e4 of the complex eigenvalues define a span formed by erotsum = (e3 + e4 ) /2 and erotdiff = (e3 − e4 ) /2i . A small displacement falls in this span will rotate arg (λ3) in the span every period of the halo orbit. Because the eigenvalue of 1 has only one eigenvector, i.e. e5 ¼ e6, the algebraic multiplicity of the monodromy matrix is two but the geometric multiplicity is one. The monodromy matrix is defective. To complete the eigenspace, the generalized eigenvectors are derived as follows, 2 el6 = Ker ⎡⎣ ( M − λ 6 I) ⎤⎦
(A.1)
Because the algebraic multiplicity is two, there are two solutions to Eq. (A.1). To have geometric meaning, the vector e˜ 6 that ^ and perpendicular to e5 is selected to lies in the span of the two e 6 ˜ be the sixth eigenvector. e6 points to nearby periodic orbits. The normalized eigenvectors e¯ i (i ¼ 1,…,6) for e1, e2, erotdiff , erotdiff , e5 and e˜ 6 form the eigenspace as follows,
E = ⎡⎣ e¯ 1 e¯ 2 e¯ 3 e¯ 4 e¯ 5 e¯ 6 ⎤⎦
(A.2)
An infinitesimal displacement of magnitude c1 along e¯ 1 will grow λ1 times in a period. Small displacements along other eigenvectors e¯ i (i¼ 2,…, 6) will not diverge as the modulus of the corresponding eigenvalues are smaller than or equal to 1. Then, as long as the components of the perturbation along each eigenvector is small,
X (X 0 + c1e¯ 1, t ) ≈ X (X 0 + EC, t )
(A.3)
where X0 is the initial state vector (manifold insertion point along the halo orbit), X is the state vector at time t, and C is the coefficient vector T C = ⎡⎣ c1 c2 c3 c4 c5 c6 ⎤⎦
(A.4)
Due to this feature, for the c1 of a given manifold, one can compute the coefficients c2,…,c6 to zero position perturbations and therefore find the Δv to approximate the effect due to sixdimensional perturbation c1e¯ 1. The constraint is expressed as,
⎡ 0(3 × 1) ⎤ ⎢ ⎥ = EC ⎣ Δv(3 × 1) ⎦
(A.5)
The matrix E is partitioned into four submatrices as,
⎡ E11 (3 × 3) E12 (3 × 3) ⎤ ⎥ E=⎢ ⎢⎣ E21 (3 × 3) E22 (3 × 3) ⎥⎦ Substituting Eq. (A.6) to Eq. (A.5) results in,
(A.6)
H. Chen et al. / Acta Astronautica 127 (2016) 464–473
⎡ c4 ⎤ ⎡ c1 ⎤ 0 = E11⎢ c2 ⎥ + E12 ⎢ c5 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ c6 ⎦ ⎣ c3⎦
References
(A.7)
and
⎡ c4 ⎤ ⎡ c1 ⎤ Δv = E21⎢ c2 ⎥ + E22 ⎢ c5 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ c6 ⎦ ⎣ c3⎦
(A.8)
Integrating Eqs. (A.7) and (A.8) to cancel [c4 c5 c6]T yields the expression for Δv; that is,
⎡ c1 ⎤ Δv = A (3 × 3) ⎢ c2 ⎥ ⎢ ⎥ ⎣ c3⎦
(A.9)
where, −1 A = (E21 − E22 E12 E11)
(A.10)
c1 is associated to the desired the manifold, and the combination of c2 and c3 determines the magnitude of Δv. Let the entry in the ith row and j-th column of A be defined as aij. Based on Eq. (A.9), the square of ΔV can be expressed as,
Δv 2 = b1c12 + b2 c22 + b3 c32 + b12 c1c2 + b23 c2 c3 + b13 c1c3
(A.11)
where, 3
bj =
3
∑ aij bkj = 2 ∑ aik aij i=1
i=1
(A.12)
At the minimum, the derivatives of Eq. (A.11) with respect to c2 and c3 should be zero; that is,
∂ΔV 2/∂c2 = b12 c1 + 2b2 c2 + b23 c3 = 0
(A.13)
∂ΔV 2/∂c3 = b13 c1 + b23 c2 + 2b3 c3 = 0
(A.14)
Then, c2 and c3 can be calculated from, −1 ⎡ c2 ⎤ ⎡ 2b2 b23 ⎤ ⎡ − b12 ⎤ ⎥ ⎢ ⎥ c1 ⎢⎣ c ⎥⎦ = ⎢ 3 ⎣ b23 2b3⎦ ⎣ − b13 ⎦
(A.15)
Let B(2×1) represent the resulting vector before c1. Then, the minimum required Δv corresponding to the manifold of c1 is expressed as,
⎡1 ⎤ Δv = A ⎢ ⎥c1 ⎣ B (2 × 1) ⎦ Meanwhile, the optimal direction to place folds is,
⎡1 ⎤ Δl v = A⎢ ⎥ ⎣ B (2 × 1) ⎦
473
(A.16)
Δv to generate mani-
(A.17)
This ΔV direction corresponds to the 6-dimensional divergent direction of the halo orbit.
[1] R.W. Farquhar, The flight of ISEE-3/ICE : origins, mission history, and a legacy, in: Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 1998. doi: 〈http://dx.doi.org/10.2514/6.1998-4464〉. [2] E.A. Belbruno, J.K. Millert, Sun-perturbed Earth-to-Moon transfers with ballistic capture, J. Guid. Control. Dyn. 16 (1993) 770–775. [3] W.S. Koon, M.W. Lo, J.E. Marsden, S.D. Ross, Low energy transfer to the Moon, Dyn. Nat. Artif. Celest. Bodies (2001) 63–73, http://dx.doi.org/10.1007/ 978-94-017-1327-6_8. [4] K.C. Howell, M. Kakoi, Transfers between the Earth–Moon and Sun–Earth systems using manifolds and transit orbits, Acta Astronaut. 59 (2006) 367–380, http://dx.doi.org/10.1016/j.actaastro.2006.02.010. [5] G. Mingotti, F. Topputo, F. Bernelli-Zazzera, Optimal low-thrust invariant manifold trajectories via attainable sets, J. Guid. Control. Dyn. 34 (2011) 1644–1656, http://dx.doi.org/10.2514/1.52493. [6] C. Martin, B.A. Conway, Optimal low-thrust trajectories using stable manifolds, in: B.A. Conway (Ed.), Spacecraft Trajectory Optimization, 1st ed.,Cambridge University Press, 2010, pp. 238–262 http://dx.doi.org/10.1017/CBO9780511778025.010. [7] D.C. Folta, M. Woodard, K. Howell, C. Patterson, W. Schlei, Applications of multi-body dynamical environments: the ARTEMIS transfer trajectory design, Acta Astronaut. 73 (2012) 237–249, http://dx.doi.org/10.1016/j. actaastro.2011.11.007. [8] M. Nakamiya, H. Yamakawa, D.J. Scheeres, M. Yoshikawa, Interplanetary transfers between halo orbits: connectivity between escape and capture trajectories, J. Guid. Control. Dyn. 33 (2010) 803–813, http://dx.doi.org/10.2514/ 1.46446. [9] G. Mingotti, F. Topputo, F. Bernelli-Zazzera, Earth–Mars transfers with ballistic escape and low-thrust capture, Celest. Mech. Dyn. Astronaut. 110 (2011) 169–188, http://dx.doi.org/10.1007/s10569-011-9343-5. [10] J.J. Guzman, D.W. Dunham, P.J. Sharer, J.W. Hunt, C.J. Ray, H.S. Shapiro, et al., STEREO mission design implementation, in: Proceedings of the 20th International Symposium on Spacecraft Flight Dynamics, Annapolis, Maryland, USA, 2007. 〈http://issfd.org/ISSFD_2007/21-4.pdf〉. [11] J. Kawaguchi, H. Yamakawa, T. Uesugi, H. Matsuo, On making use of lunar and solar gravity assists in LUNAR-A, PLANET-B missions, Acta Astronaut. 35 (1995) 633–642, http://dx.doi.org/10.1016/0094-5765(95)00013-P. [12] S. Campagnola, C. Corral, R. Jehn, Design of lunar gravity assist for the bepicolombo mission to mercury, Adv. Astronaut. Sci. 119 (2004) 427–442. [13] D. Landau, T.P. Mcelrath, D. Grebow, N.J. Strange, Efficient lunar gravity assists for solar electric propulsion missions, Adv. Astronaut. Sci. 143 (2012) 917–934. [14] T.P. McElrath, G. Lantoine, D. Landau, D. Grebow, N. Strange, R. Wilson, et al., Using gravity assists in the Earth–Moon system as a gateway to the solar system, in: Proceedings of the Global Space Exploration Conference, Washington, D.C., 2012. [15] H. Chen, Y. Kawakatsu, T. Hanada, Earth Escape from a Sun-Earth halo orbit using unstable manifold and lunar swingbys, Trans. Jpn. Soc. Aeronaut. Space Sci., 59, 2016. [16] Y. Kawakatsu, Destiny mission overview: a small satellite mission for deep space exploration technology demonstration, Adv. Astronaut. Sci. 146 (2012) 727–739. [17] V. Szebehely, Theory of Orbits: The Restricted Problem of Three Bodies, Academic Press, New York, San Francisco, London 1967, pp. 134–139. [18] W.S. Koon, M.W. Lo, J.E. Marsden, S.D. Ross, Dynamical systems, Three-Body Probl. Sp. Mission Des. (2011) 174–175. [19] K.C. Howell, Three-dimensional, periodic, “Halo” orbits, Celest. Mech. (1984). [20] G. Lantoine, T.P. Mcelrath, Families of solar-perturbed Moon-to-Moon transfers, in: Proceedings of theAAS/AIAA Space Flight Mechanics Meeting, 2014. [21] R.R. Bate, D.D. Mueller, J.E. White, Fundamentals of astrodynamics, Dover (1971). [22] N. Arora, R.P. Russell, N.J. Strange, Partial Derivatives of the Lambert Problem, in: Proceedings of the AIAA/AAS Astrodyne Specialist Conference, American Institute of Aeronautics and Astronautics, 2014. 〈http://dx.doi.org/10.2514/6. 2014-4427〉. [23] H. Chen, NEW Results of Orbital Dynamics and the Application to Orbit Prediction and Mission Design, Kyushu University 2015, pp. 148–153. [24] F. Topputo, R. Zhang, F. Bernelli-Zazzera, J. Luo, Numerical approximation of invariant manifolds in the restricted three-body problem, in: Proceedings of the International Astronautical Congress, Beijing, China, 2013. [25] R. Keys, Cubic convolution interpolation for digital image processing, Acoust. Speech Signal Process. IEEE Trans 29 (1981) 1153–1160, http://dx.doi.org/ 10.1109/TASSP.1981.1163711.