- Email: [email protected]
Deep space exploration orbit design departing from circumlunar orbit of lunar base
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.1 (1-13)
Aerospace Science and Technology ••• (••••) ••••••
Contents lists available at ScienceDirect
Aerospace Science and Technology www.elsevier.com/locate/aescte
Deep space exploration orbit design departing from circumlunar orbit of lunar base Shi Qiu a , Xibin Cao a , Feng Wang a,∗ , Chengfei Yue a , Zhe Zhang b a b
Research Center of Satellite Technology, Harbin Institute of Technology, Harbin 150001, China Lunar Exploration and Space Engineering Center, China
a r t i c l e
i n f o
Article history: Received 15 May 2019 Received in revised form 11 October 2019 Accepted 18 October 2019 Available online xxxx Keywords: Circumlunar orbit Lunar base Fixed excess velocity vector Analytical trajectory design Genetic algorithm
a b s t r a c t A design method for a new kind of deep space mission which departs from a circumlunar orbit of the lunar base is presented in this paper. The difference between this kind of deep space mission and the traditional one is well analyzed. The method of patched conics is employed to divide the mission trajectory into four independent parts. The independent selenocentric segment is designed using a Lambert Problem based method, while the geocentric trajectory from a certain point in threedimensional space to a fixed excess velocity vector is determined analytically. Then the combination of these segments is optimized via genetic algorithm by minimizing the total velocity increasement. Finally, typically numerical simulations are demonstrated to verify the effectiveness of the proposed method. © 2019 Elsevier Masson SAS. All rights reserved.
1. Introduction Deep space exploration has significant importance in exploring the distant regions of outer space. In deep space exploration, orbit design and the orbital optimization [1,2] in various scenarios are usually the first and one of the most important steps before the spacecraft functions well under virous control strategy [3–5]. For traditional deep space missions departing from a low circular geocentric parking orbit, the orbit design problem and optimal maneuver has been studied intensively [6], and feasible impulsive transfers to some fixed hyperbolic excess velocity vector have been obtained [7,8]. The obtained impulsive solutions can be used as the initial values to determine finite-burn or continuous transfers. Among these impulsive solutions, the three-impulse time-free initial guess method reflects the behavior of the optima, resulting in better overall convergence, while the three-impulse time-fixed one and optimization method allows for quick cost assessments for a variety of escape geometries as a function of the total trip time [9,10]. However, these impulsive based methods generate orbital error in the practical application, which should be adjusted by trajectory correction maneuver using some strategies, for instance, the B-plane targeting method [11].
*
Corresponding author. E-mail address: [email protected] (F. Wang).
https://doi.org/10.1016/j.ast.2019.105505 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.
For the deep space exploration, lunar exploration and the exploration from lunar base has attracted a great deal of attention. Many countries have established different exploration missions, such as the Apollo Missions, Chang’e Missions. In these missions, the earth-to-moon transfer orbit and the return orbit design is the preliminary task. Many scholars have studied the earth-to-moon transfer orbits and the return orbits, and the design results have been verified by practical missions [12–16]. For some unexpected conditions such as a brake failure at the perigee of the moon, the free return orbits offers an optional way of coming back [17–21], and the application in Appollo-13 was called “a successful failure” [22]. Based on these remarkable results, regenerative and sustainable hybrid life support systems for a lunar base has even been designed [23]. All these improvements in lunar exploration pave the way of deep space exploration departs from lunar base. In this paper, a design method for a new kind of deep space mission which departs from a circumlunar orbit of lunar base to another planet is presented. Firstly, the difference to a traditional deep space mission starts from a low geocentric orbit is analyzed in section 2. Following that in section 3, the selenocentric trajectory is calculated on the basis of Lambert problem. The formation of initial position and velocity vectors of the geocentric trajectory are also obtained. In section 4, a new orbit design problem that departs from a certain point to a known excess velocity is presented. In order to solve this problem, two auxiliary frame system were established, thus analytic solutions are derived from a sixth-degree
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.2 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
2
equation. Next, Newton-Raphson iteration algorithm is applied to make the whole trajectory time-continuous and the genetic algorithm is used to find the best solution in section 5. Examples of deep space missions from a circumlunar orbit are demonstrated to validate the effectiveness in section 6. Finally, section 7 provides concluding remarks of this paper. 2. Problem description Launch opportunity determination is the first step of a deep space orbit design. When the ranges of launch time and flight time are settled, the best launch opportunity with the smallest total velocity increasement during the interplanetary transferring can be found using pork-chop method [24]. Once the launch opportunity is settled, the excess velocity vector denoted by v ∞ and the corresponding time instant t ∞ can be determined as well. A traditional deep space exploration mission usually starts from a low geocentric parking orbit. In these missions, when the excess velocity, the height of the acceleration point and orbital inclination are determined, the B-plane method can be applied to calculate the position and the velocity of the acceleration point. Then the geocentric hyperbolic orbit from the acceleration point to the excess velocity is determined. For the mission starting from a circumlunar orbit, the design method is quite different. The patched conics method is employed to divide the mission trajectory into four parts: the departure trajectory relative to the moon; the hyperbolic orbit relative to the earth, the cruise ellipse relative to the sun, and the hyperbolic arrival trajectory relative to the target planet. Firstly, the spacecraft is mainly influenced by the gravity of the moon, and it has to get rid of the lunar gravity in order to get into the sphere of influence of the earth, as shown in Fig. 1. The connection point of the selenocentric segment and the geocentric segment is the exit point through which the spacecraft flies out of the sphere of influence of the moon. This exit point is naturally the beginning of the geocentric hyperbolic orbit. However, this exit point can be arbitrary on the lunar sphere of influence, which leads to a geocentric orbit design problem starting from a certain point rather than a parking orbit to an excess velocity. Once this geocentric hyperbolic orbit with the excess velocity vector v ∞ is determined, the rest parts of this mission is similar to a traditional mission. Thus, the last two parts of the orbit in this mission will not be discussed in detail in this article. So, the purpose of this paper is to find an optimal orbit with a fixed excess velocity optimized over the exit point on the influence sphere of the moon as illustrated in Fig. 1. 3. Selenocentric segment design method In order to combine the selenocentric segment with the geocentric segment easily, the moon orbital plane frame coordinate O L -xyz is established [25]. The origin O L is the center of the moon, and the x- y plane is the moon’s orbital plane at the time instant t B when the spacecraft arrives at the boundary of the lunar sphere of influence. The x axis directs from the moon towards the earth while the z axis is normal to the moon’s orbital angular momentum vector. In addition, the latitude of the lunar orbit is set to zero while x-z plane is the meridian plane as is shown in Fig. 2. The point where the spacecraft flies out of the sphere of the lunar influence is denoted by B. Assume that the latitude of B is ϕ B while the longitude is λ B . Thus, the position vector of exit B can be expressed as
⎡
⎤
cos ϕ B · cos λ B L r BL = ⎣ cos ϕ B · sin λ B ⎦ · ρsoi sin ϕ B
(1)
Fig. 1. Deep space mission from a selenocentric orbit.
Fig. 2. Moon orbital plane frame coordinate. L where ρsoi is the radius of sphere of influence of the moon, and the specific value is 66200 kilometers. In this mission all six initial orbital elements are set, then the position information of the spacecraft on the parking orbit can be acquired easily at any given time with the help of high-precision gravitational model. In addition, the position of exit B can be obtained by equation (1). Once the trip time from the parking orbit to the exit is acquired, by solving the Lambert problem about T L (which will be discussed in section 5), the orbit of selenocentric segment is determined. According to the definition of the moon orbital plane frame coordinate O L -xyz, it is easy to calculate its orthogonal unit basis vectors with respect to the earth-J2000 system as follows
L i L = r E , |r EL |
k L =
r EL × v LE , |r EL × v LE |
j L = k L × i L
(2)
where, r EL and v LE are the moon’s position and the velocity relative to the earth at time instant t B when the spacecraft arrives at the exit of lunar sphere of influence, and they can be both obtained by searching high-precision ephemeris, DE421 for example. Thus, the transformation from moon orbital plane frame coordinate O L -xyz to the earth-J2000 system yields:
R LE = [i L , j L , k L ]
(3)
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.3 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
3
Fig. 4. Hyperbolic orbit description frame.
√ sin |β| =
Fig. 3. Hyperbolic orbit description frame.
Therefore, the end position of the selenocentric segment which is the initial position of the geocentric segment relative to the earth can be expressed as follows
r BE = R LE · r BL + r LE
(4)
4. Geocentric segment design method Till now, the initial position r BE and hyperbolic excess velocity v ∞ with respect to the earth-J2000 frame are obtained. Then the geocentric orbit design problem from a certain position r BE to an excess velocity vector v ∞ will be discussed in this section. As shown in Fig. 3, the earth lies at one of the focus of hyperbolic orbit. In order to describe the orbit easily, we established a new reference frame H-xyz. The x- y plane of the H-xyz frame is the plane of the orbit, and the x axis is directed along the apse line, and the center H is |a| + r p away from the earth. The z axis is normal to the plane of the orbit in the direction of the orbital angular momentum vector. In this new frame, the orbit can be described in hyperbolic equation as
x2
|a|2
y2
−
|b|2
=1
(5)
where a is the semimajor axis while b is the semiminor axis that can be found in terms of a and the eccentricity e of the hyperbolic orbit as
b = a e2 − 1
(6)
According to [26], the specific energy of a hyperbolic orbit is
ε=
μE
(7)
2|a|
with μ E the gravitational parameter of the earth. Substituting ε , v ∞ into to living force equation yields
|v ∞ |2 2 where
|a| =
−
μE μE = E 2|a| ρsoi
(8)
E ρsoi is the radius of sphere of influence of the earth. Thus,
μE |v ∞ |2 − 2ρμE E
(9)
soi
The angle β in Fig. 3 is the angle between asymptote and apse, which satisfies
e2 − 1
cos |β| =
,
e
1
(10)
e
To get the transformation from the earth-J2000 frame to H-xyz system, an auxiliary frame system M-xyz is established. The x axis of M-xyz is directed along the asymptote of the hyperbola and the z axis is normal to the plane of the orbit in the direction of the orbital angular momentum vector, as is shown in Fig. 4. Considering that the excess velocity is parallel to the asymptote when they approach infinity, the orthogonal unit basis vectors for this auxiliary frame can be expressed as follows
i M = v ∞ , |v ∞ |
k M =
r0 × v ∞ , |r0 × v ∞ |
j M = k M × i M
(11)
Then the transformation from the earth-J2000 frame to the Mxyz system can be obtained as:
T T T RM J = i M , j M , kM
(12)
Then the initial position vector r BE in the M-xyz frame system is
r M = R MJ · r BE
(13)
By now, it is easy to obtain the transformation from M-xyz to H-xyz by rotating a angle β about the z axis. The sign of angle β is positive when the rotation is clockwise. Otherwise, the angle β is negative. Thus, transformation from M-xyz to H-xyz is
⎡
H RM
cos β = ⎣ − sin β 0
⎤
sin β cos β 0
0 0⎦ 1
⎧⎡ e 2 −1 1 ⎪ ⎪ e e ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ e 2 −1 1 ⎪ ⎣− e ⎪ ⎪ e ⎪ ⎪ ⎨ 0 0 = ⎡ 2 ⎪ 1 ⎪ ⎪ − e e −1 ⎪ e ⎪ ⎢ ⎪ ⎪ ⎢ 2 ⎪ 1 ⎪ ⎣ e −1 ⎪ e e ⎪ ⎩
0
0
⎤ 0
⎥
⎥ 0⎦ β >0 1 0
⎤
(14)
⎥
⎥ 0⎦ β <0 1
Therefore, the transformation from the earth-J2000 to H-xyz can be calculated by (12) and (14): H R EH = R M · RM E
(15)
Furthermore, the origin of H-xyz system is |a| + r p behind the origin of the earth-J2000 system. The initial position r BE relevant the H-xyz system can be expressed as follows
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.4 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
4
⎡ H r H = R M · RM r BE + ⎣ E ·
|a| + r p 0 0
⎤
⎡
|a|e
⎤
H ⎦ = RM · r M + ⎣ 0 ⎦
(16)
0
where r p is the perigee radius of the hyperbolic orbit that can be found in terms of a and e as
r p = |a|(e − 1) Substituting
(17) H RM
into to equation (16) yields
⎧⎡ ⎤ e 2 −1 1 M ⎪ ⎪ rx + e r M + |a|e y ⎪ e ⎪ ⎥ ⎢ ⎪ ⎪ ⎥ β >0 ⎢ ⎪ e 2 −1 M 1 M ⎪ ⎦ ⎣ − r + r ⎪ x e e y ⎪ ⎪ ⎪ ⎨ 0 r H = ⎡ ⎤ 2 ⎪ e −1 M 1 M ⎪ ⎪ r − r + | a | e y ⎪ e ⎪ ⎢ e x ⎥ ⎪ ⎪ ⎢ ⎥ β <0 ⎪ e 2 −1 M ⎪ ⎣ ⎦ r x + 1e r M ⎪ y e ⎪ ⎪ ⎩
(18)
0
Fig. 5. Example of multiple solutions and restricted area.
where r xM and r M r M . The hypery are the first two components of bolic orbit lies in the x- y plane of the H-xyz system, so the third component of r M is zero. Taking the case β > 0 as example, substituting r H into to the hyperbolic equation, we can obtain
( 1e r xM +
e 2 −1 M ry e a2
Denoting x =
+ |a|e )2
−
(−
√
e 2 − 1 yields
2 r xM + r M y x + |a| x − 1
2
e 2 −1 M 2 r x + 1e r M y ) e 2 2 a (e − 1)
x2 − −r xM x + r M y
2
=1
(19)
= |a|x2 x2 − 1 (20)
The above equation can further be written in the following formation. 6
5
4
3
2
c6 x + c5 x + c4 x + c3 x + c2 x + c1 x + c0 = 0
(21)
Therefore, a sixth-degree closed form equation with only one variable x is obtained. All the coefficients of this equation are expressed as follows:
⎧ ⎪ c = a2 ⎪ ⎪ 6 ⎪ ⎪ ⎪ c 5 = 2|a|r M ⎪ y ⎪ ⎪ ⎪ M 2 ⎪ ⎪ ⎪ c = ry + 2|a|r xM + a2 ⎪ ⎨ 4 M c 3 = 2r xM r M y + 2|a|r y ⎪ ⎪ ⎪ c = 2|a|r M ⎪ ⎪ 2 x ⎪ ⎪ ⎪ M M ⎪ ⎪ c 1 = 2r x r y ⎪ ⎪ ⎪ ⎪ ⎩ c = −r M 2 0
r x2
−
a2
r 2y a2 (e 2 − 1)
H r Bx
a2
x−
y
(23)
(24)
r BHy a2 (e 2 − 1)
y=1
(25)
Based on the energy conservation principle, the magnitude of velocity at exit B can be obtained from the following equation.
1 2
v 2B −
μE |r B |
=
1 2
v 2∞ −
μE |r∞ |
(26)
Then the velocity | v B | can be obtained:
|v B | =
(22)
=1
According to tangent formula of hyperbolic, the direction of velocity at exit B in geocentric segment must be parallel to the following line:
Similarly, when β is negative, we can also acquire a closed form equation of sixth degree. The coefficients are
⎧ ⎪ c 6 = a2 ⎪ ⎪ ⎪ ⎪ ⎪ c 5 = −2|a|r M ⎪ y ⎪ ⎪ ⎪ M 2 ⎪ ⎪ ⎪ c = ry + 2|a|r xM + a2 ⎪ ⎨ 4 M c 3 = −2r xM r M y − 2|a|r y ⎪ ⎪ ⎪ ⎪ c 2 = 2|a|r xM ⎪ ⎪ ⎪ ⎪ M M ⎪ ⎪ ⎪ c 1 = −2r x r y ⎪ ⎪ ⎪ ⎩ c = −r M 2 0 y
By solving the equation, the eccentricity of the hyperbola can be obtained, so that the hyperbolic orbit equation in the H-xyz system yields
v 2∞ +
2μ E
|r B |
−
2μ E
|r∞ |
(27)
It is worth noticing that the direction of the orbital angular momentum is relevant to flight direction of the orbit. Therefore, there may be more than one orbits, one clockwise hyperbolic orbit and an anticlockwise one, departing from a certain position to a given excess velocity vector satisfying the boundary conditions. What’s more, for each certain excess velocity v ∞ , there is a limited area where only one hyperbolic orbit can be found satisfying given conditions. Take Fig. 5 as an example, the left boundary of the shadow region is parallel to the excess velocity vector while the right one is symmetric with the left one about x axis. Therefore, the shadow region will be feasible region of the anticlockwise hyperbolic obits with the earth as their focus and the exclusion region of the clockwise orbits. However, solving the equation (21) may result in some infeasible solutions because the excess velocity vector only represents the direction not the position when the spacecraft flies out of the influence sphere of the earth. For example, the hyperbolic orbit 3 in Fig. 5 is an infeasible solution passing through the point P and with a same velocity magnitude of v ∞ . In order to reject these infeasible solutions, we can calculate the angle γ between the excess
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.5 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
5
In order to make the entire trajectory time continuous, the time t B has to be chosen consistently such that: t∞ = t∞
(31)
Then t B can be chosen via Newton-Raphson iteration according to [27]. Up to now, the geocentric trajectory along with the heliocentric one has been determined. According to the analysis of section 1, once the trip time in the lunar sphere of influence is settled, the selenocentric trajectory can be acquired by solving the Lambert problem. Thus, the departure velocity and the velocity of the arrival at the exit are both obtained. Therefore, the total velocity increasement yields
− v = |v d − v 0 | + v + B − vB
Fig. 6. Example of virtual branch of hyperbolic orbit.
velocity vector v ∞ and the position vector r EA , in which E represents the earth while A represents the exit point of the earth’s influence sphere.
γ = arc cos
v ∞ · r EA
(28)
If γ is an acute angle, this solution is a feasible one. Otherwise when γ is an obtuse angle, the solution will be infeasible. Considering the symmetricity of the hyperbola, the real trajectory is always accompanied with its reflection. Mathematically, if we replace x and y with −x and − y, we will obtain the same equation:
(−x)2 (− y )2 x2 y2 − = − =1 |a|2 |b|2 |a|2 |b|2
(29)
Thus, the solution to the equation (21) may also result in a virtual trajectory as is shown in Fig. 6, which forms the infeasible hyperbolic orbit 3 in Fig. 5. 5. Optimization of the orbit When the pork-chop method is applied to finding the best opportunities to travel from the earth to another planet, the time t ∞ of the excess velocity v ∞ is determined. Meanwhile, the solution to the geocentric segment is only depend on the initial position vector r BE and the excess velocity v ∞ , which means the trip time of the geocentric segment denoted by T E is fixed when r BE and v ∞ are chosen. For the initial position vector r BE , it is determined by the latitude and longitude of the exit of lunar sphere of influence as well as the moon’s real-time position vector r LE in the sphere of influence of the earth. Hence, if the time t B when the spacecraft arrives at the exit of lunar sphere of influence is determined, the real time that the spacecraft flies out of the sphere of influence of the earth is settled as well. t∞ = tB + T E
where, v d represents the velocity of departure while v 0 is the velocity of the parking orbit. v + B is the initial velocity of the geocentric trajectory while v − B is the ending velocity of the selenocentric trajectory. Since the total velocity increasement is not an analytical solution to the trip time of the selenocentric segment T L , an intelligent optimization algorithm can be employed to minimize the total velocity increasement with respect to T L . Then the optimization index, the total velocity increasement, can be expressed by a function of T L as:
− J = v = | v d − v 0 | + v + B − v B = f (T L )
|v ∞ | · |r EA |
(32)
(30)
where, t ∞ is the real time that the spacecraft flies out of the sphere of influence of the earth.
(33)
To obtain the optimal solution of equation (33) the genetic algorithm, owing to its significant advantage of the capability of initial-condition-independent global optimization and the strong robustness, has been adopted. The procedure of solving this optimization problem is listed as follows and Fig. 7. Step 1. Using pork-chop method to find the best opportunity departing from the earth to another planet, i.e., to find the excess velocity vector v ∞ and its time t ∞ ; Step 2. By choosing the latitude λ B and longitude ϕ B of the exit of lunar sphere of influence, as well as the time t B , to get the initial position relative to the earth r BE ; Step 3. Using the results of v ∞ and r BE and the method introduced in section 3 to get geocentric segment; Step 4. Calculate the trip time from the initial position to the excess velocity vector v ∞ ; and verify whether it Step 5. Calculate the real arrived time t ∞ satisfies the time error condition or not; Step 6. Repeat Step 2 to Step 5 until the real arrived time t ∞ satisfies the time error condition; Step 7. Calculate the optimal total velocity increasement. 6. Simulation cases In order to validate the effectiveness of the design method, three simulation cases with different target planets are presented in this section. For the first two examples, Mars exploration missions are considered while the third one considers the Venus exploration mission. Thus, the effectiveness of proposed method with different launching time and different target planets are both verified through the listed scenarios. Example 1 Assume a deep space exploration mission to the mars which will be carried out from a circumlunar orbit during 2020. The initial orbital elements relevant to the moon − J 2000 frame of a circumlunar spacecraft are shown in Table 1.
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.6 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
6
Fig. 7. Flow chart of the mission orbit design. Table 1 Initial orbital elements of a circumlunar orbit. Time (UTCG)
Semimajor axis (km)
Eccentricity
Inclination (deg)
RAAN (deg)
Arg of perigee (deg)
True anomaly (deg)
1 Jul 2020 00:00:00.000
8943.996897
0.693000
25.400
275.800
182.400
341.200
Fig. 9. Three-dimensional figure of total velocity increasement in 2020. Fig. 8. Pork-chop figure of total velocity increasement in 2020.
According to the analysis in this article, the excess velocity vector v ∞ and its time t ∞ have to be calculated at first. A pork-chop figure about total velocity increasement is shown in Fig. 8 (Fig. 9). It is easy to find the best opportunity using an optimal method in this pork-chop figure, marked with a star symbol “∗”. The
departure time and the excess velocity of the best opportunity yield
t ∞ = 2459054.7534 (Julian day) v ∞ = [3.2298; 0.7070; 1.6406] (km/s)
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.7 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
7
Table 2 Coefficients of equation (21) when β is positive. Coefficients
c 0 (e+11)
c 1 (e+10)
c 2 (e+9)
c 3 (e+9)
c 4 (e+11)
c 5 (e+10)
c 6 (e+8)
Values
-1.1620
1.1752
-1.0769
-9.5428
1.1610
-2.1295
9.7565
Table 3 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
11.8796
9.8612
1.0521
0.0+1.0i
0.0-1.0i
-0.9663
Table 4 Coefficients of equation (21) when β is negative. Coefficients
c 0 (e+11)
c 1 (e+10)
c 2 (e+9)
c 3 (e+9)
c 4 (e+11)
c 5 (e+10)
c 6 (e+8)
Values
-1.1620
-1.1752
-1.0769
9.5428
1.1610
2.1295
9.7565
Table 5 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
-11.8796
-9.8612
-1.0521
0.0+1.0i
0.0-1.0i
0.9663
Table 6 Orbital elements of feasible solutions of geocentric segment. Conditions
Julian day of exit B
Eccentricity of geocentric segment
r B in earth-J2000 (km)
v B in earth-J2000 (km/s)
β >0
2459052.0256 2459052.0369
11.7708 9.7677
[−169592.6; 257181.3; 135567.6] [−170361.7; 256529.5; 135359.7]
[3.2529; 1.0112; 1.8777] [3.229; 1.057; 1.8995]
β <0
2459051.1055
1.3070
[−101775.2; 303652.0; 149094.7]
[0.8726; −3.3514; −1.7391]
To a common deep space mission, it usually takes about 3 days to fly out of the sphere of influence of the earth. Thus, the initial guessed value of t B can be assumed to be 3 days earlier than t ∞ . Meanwhile, in this example case, the latitude and longitude of the exit of lunar sphere of influence are both set to zero. When β is positive, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 2. It is easy to obtain the solutions to the six-degree equation which are listed in Table 3. It is apparently that x4 , x5 and x6 are infeasible solutions because the eccentricity have to be a real value and must be greater than 1. Therefore, potential feasible solutions must be among x1 , x2 and x3 . As for the solution x1 = 11.8796, initial trip time is 2.6581 days, not equal to 3 days. Using Newton-Raphson iteration algorithm, we can obtain convergent value, x1 = 11.7282, under which condition, the trip time of geocentric segment is 2.7278 days. Similarly, for the solution x2 = 9.8612, the convergent value is x2 = 9.7164 and the geocentric trip time is 2.7165 days. However, for the solution x3 = 1.0521, the excess velocity vector v ∞ relative to H-xyz frame is [2.5428; −2.6753; 0.0000] while the velocity when the spacecraft flies out of the influence sphere of the earth is [−2.5428; −2.6753; 0.0000]. That means this solution lies on the other branch of the hyperbolic orbit which does not exit indeed. When β is negative, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 4. It is easy to obtain the solutions to the six-degree equation which are listed in Table 5. It is obviously that only x6 may be a feasible solution. The convergent value is x1 = 0.8415 and the geocentric trip time is 3.6479 days.
Fig. 10. Convergence curve of the first case when β is positive.
The orbital elements of feasible solutions are listed in Table 6. Thus, the geocentric trajectory has been obtained. The next step is to find a best selenocentric solution to minimize the total velocity increasement the least according to the optimal method in section 5. Genetic algorithm is applied in this case. The convergence curves under different conditions are shown in Figs. 10 to 12. It is obvious that genetic algorithm is suitable in this optimal problem for its fast convergence. The best solutions are illustrated in Table 7.
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.8 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
8
Table 7 Orbital elements of selenocentric segment. Conditions
Selenocentric trip time
Acceleration position in moon-J2000 (km)
Acceleration velocity in moon-J2000 (km/s)
β >0
1.9315 days 2.8296 days
[−4636.9; −8371.5; −2592.2] [−4748.4; −7881.5; −2621.3]
[0.2788; −0.8237; −0.3637] [0.1750; −0.8481; −0.3698]
β <0
0.5507 days
[5639.9; −12133.7; 2082.1]
[−0.7089; 1.5064; −0.1930]
Fig. 13. Pork-chop figure of total velocity increasement in 2022.
Fig. 11. Convergence curve of the second case when β is positive.
Fig. 12. Convergence curve when β is negative.
Example 2 Assume another deep space exploration mission to the mars which will be carried out from a circumlunar orbit during 2021. The initial orbital elements relevant to the moon − J 2000 frame of a circumlunar spacecraft are shown in Table 8.
Fig. 14. Three-dimensional figure of total velocity increasement in 2022.
The pork-chop figure about total energy consumption in this case is shown in Fig. 13 (Fig. 14).
Table 8 Initial orbital elements of a circumlunar orbit. Time (UTCG)
Semimajor axis (km)
Eccentricity
Inclination (deg)
RAAN (deg)
Arg of perigee (deg)
True anomaly (deg)
1 Sep 2022 00:00:00.000
8943.996897
0.693000
25.400
275.800
182.400
341.200
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.9 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
9
Table 9 Coefficients of equation (21) when β is positive. Coefficients
c 0 (e+10)
c 1 (e+11)
c 2 (e+10)
c 3 (e+11)
c 4 (e+10)
c 5 (e+9)
c 6 (e+8)
Values
-3.8284
-1.4887
1.4648
-1.5640
5.3303
-7.5338
3.7064
Table 10 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
10.3998
5.0816+4.0197i
5.0816-4.0197i
0.0+1.0i
0.0-1.0i
-0.2366
Table 11 Coefficients of equation (21) when β is negative. Coefficients
c 0 (e+10)
c 1 (e+11)
c 2 (e+10)
c 3 (e+11)
c 4 (e+10)
c 5 (e+9)
c 6 (e+8)
Values
-3.8284
1.4887
1.4648
1.5640
5.3303
7.5338
3.7064
Table 12 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
-10.3998
5.0816+4.0197i
5.0816-4.0197i
0.0+1.0i
0.0-1.0i
0.2366
Table 13 Orbital elements of feasible solutions of geocentric segment. Conditions
Julian day of exit B
Eccentricity of geocentric segment
r B in earth-J2000 (km)
v B in earth-J2000 (km/s)
β >0 β <0
2459840.6552 2459838.8930
12.2782 1.0274
[83290.5; 359056.7; 182649.6] [225503.9; 329137.8; 157305.9]
[2.4096; 2.2604; 3.4244] [−2.5021; −3.6221; −1.7841]
It is easy to find the best opportunity using an optimal method in this pork-chop figure, marked with a star symbol “∗”. The time and the excess velocity of the best opportunity yield
t ∞ = 2459530.0293 v ∞ = [0.3589; −4.1368; 0.7600] Also, in this example, the initial guessed value of t B is assumed to be 3 days earlier than t ∞ and the latitude and longitude of the exit of lunar sphere of influence are both set to zero. When β is positive, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 9. It is easy to obtain the solutions to the six-degree equation which are listed in Table 10. Obviously, only x1 can be a potential feasible solution. The convergent value is x1 = 12.2372 and the geocentric trip time is 1.3679 days. When β is negative, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 11. It is easy to obtain the solutions to the six-degree equation which are listed in Table 12. It is obviously that only x6 may be a feasible solution. The convergent value is x6 = 0.2357 and the geocentric trip time is 3.1301 days. The orbital elements of feasible solutions are listed in Table 13. Thus, the geocentric trajectory has been obtained. The next step is to find a best selenocentric solution to make the total velocity increasement the least according to the optimal method in section 5. Genetic algorithm is applied in this case. The convergence curves under different conditions are shown in Fig. 15 and Fig. 16.
Fig. 15. Convergence curve when β is positive.
It is obvious that genetic algorithm is suitable in this optimal problem for its fast convergence speed. The best solutions are displayed in Table 14. Nevertheless, it is worth noticing that, in the case β is negative, the spacecraft will fly by the earth rather than directly fly out of the sphere of influence of the earth when β is positive. Meanwhile, the perigee radius of the hyperbolic orbit in the negative case is smaller than the radius of the earth. That means, the spacecraft will intersect into the earth, as shown in Fig. 17. Therefore, the trajectory in the negative case cannot be applied into engineering.
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.10 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
10
Example 3 Assume another deep space exploration mission to the Venus which will be carried out from a circumlunar orbit during 2023. The initial orbital elements relevant to the moon − J 2000 frame of a circumlunar spacecraft are shown in Table 15. The pork-chop figure about total energy consumption in this case is shown in Fig. 18 (Fig. 19). It is easy to find the best opportunity using an optimal method in this pork-chop figure, marked with a star symbol “∗”. The time and the excess velocity of the best opportunity yield
t ∞ = 2459530.0291 v ∞ = [−1.8218; 2.8860; −0.6494]
Fig. 16. Convergence curve when β is negative.
Also, in this example, the initial guessed value of t B is assumed to be 3 days earlier than t ∞ and the latitude and longitude of the exit of lunar sphere of influence are both set to zero. When β is positive, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 16. It is easy to obtain the solutions to the six-degree equation which are listed in Table 17.
Table 14 Orbital elements of selenocentric segment. Conditions
Selenocentric trip time
Acceleration position in moon-J2000 (km)
Acceleration velocity in moon-J2000 (km/s)
β >0 β <0
3.1925 days 0.5507 days
[6215.8; −3382.3; 2774.1] [6186.2; −3282.0; 2764.9]
[−0.5033; 0.9420; 0.0928] [−1.2829; 1.1034; −0.4854]
Fig. 17. Trajectory in earth-J2000 frame when β is negative.
Table 15 Initial orbital elements of a circumlunar orbit. Time (UTCG)
Semimajor axis (km)
Eccentricity
Inclination (deg)
RAAN (deg)
Arg of perigee (deg)
True anomaly (deg)
25 May 2023 00:00:00.000
2037.400338
0
45.398
32.335
190.077
43.793
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.11 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
11
It is apparently that x3 , x5 and x6 are infeasible solutions because the eccentricity have to be a real value and must be greater than 1. Therefore, potential feasible solutions must be among x1 , x2 and x4 . As for the solution x1 = 10.1589, initial trip time is 3.2148 days, not equal to 3 days. Using Newton-Raphson iteration algorithm, we can obtain convergent value, x1 = 10.2052, under which condition, the trip time of geocentric segment is 3.1749 days. Similarly, for the solution x4 = 0.7054, the convergent value is x4 = 0.8594 and the geocentric trip time is 3.8180 days. However, for the solution x2 = 8.1155, the excess velocity vector v ∞ relative to H-xyz frame is [0.4249; −3.4481; 0.0000] while the velocity when the spacecraft flies out of the influence sphere of the earth is [0.4249; 3.4481; 0.0000]. That means this solution lies on the other branch of the hyperbolic orbit which does not exit indeed. When β is negative, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 18. It is easy to obtain the solutions to the six-degree equation which are listed in Table 19. It is obviously that only x3 may be a feasible solution. The convergent value is x3 = 1.1401 and the geocentric trip time is 3.8257 days. The orbital elements of feasible solutions are listed in Table 20. Thus, the geocentric trajectory has been obtained. The next step is to find a best selenocentric solution to make the total velocity increasement the least according to the optimal method in section 5. Genetic algorithm is applied in this case. The convergence curves under different conditions are shown in Figs. 20 to 22. It is obvious that genetic algorithm is suitable in this optimal problem for its fast convergence speed. The best solutions are displayed in Table 21.
Fig. 18. Pork-chop figure of total velocity increasement in 2023.
7. Conclusion In this paper, we introduced a new type of deep space mission departing from a circumlunar orbit. Since the initial central body is the moon, not the earth, the selenocentric trajectory is a new part of the whole mission compared with the traditional one. The mission trajectory is divided into four independent parts using patched conics method. The method based on Lambert problem
Fig. 19. Three-dimensional figure of total velocity increasement in 2023. Table 16 Coefficients of equation (21) when β is positive. Coefficients
c 0 (e+10)
c 1 (e+10)
c 2 (e+9)
c 3 (e+10)
c 4 (e+10)
c 5 (e+10)
c 6 (e+9)
Values
-9.8415
8.7773
-9.9506
6.5459
8.9729
-2.2314
1.2648
Table 17 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
10.1589
8.1155
-1.3379
0.7054
+1.0i
-1.0i
Table 18 Coefficients of equation (21) when β is negative. Coefficients
c 0 (e+10)
c 1 (e+10)
c 2 (e+9)
c 3 (e+10)
c 4 (e+10)
c 5 (e+10)
c 6 (e+9)
Values
-9.8415
-8.7773
-9.9506
-6.5459
8.9729
2.2314
1.2648
Table 19 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
-10.1589
-8.1155
1.3379
-0.7054
+1.0i
-1.0i
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.12 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
12
Table 20 Orbital elements of feasible solutions of geocentric segment. Conditions
Julian day of exit B
Eccentricity of geocentric segment
r B in earth-J2000 (km)
v B in earth-J2000 (km/s)
β >0
2460098.3879 2460097.7448
10.2541 1.3185
[−147868.0; −277401.4; −134052.0] [−194745.1; −248150.1; −116028.6]
[−2.2208; 2.7903; −0.9073] [2.3835; 2.4533; 1.3712]
β <0
2460097.7371
1.5165
[−195271.9; −247764.1; −115795.8]
[2.4056; 2.4296; 1.3749]
Table 21 Orbital elements of selenocentric segment. Conditions
Selenocentric trip time
Acceleration position in moon-J2000 (km)
Acceleration velocity in moon-J2000 (km/s)
β >0
1.6314 days 0.5991 days
[−615.6; −1627.2; −1060.2] [262.1; −1461.9; −1394.6]
[1.0030; −0.2403; −0.1578] [1.5493; 0.2189; 0.0565]
β <0
0.5913 days
[274.4; −1457.2; −1397.2]
[1.5619; 0.2283; 0.0609]
Fig. 20. Convergence curve of the first case when β is positive.
Fig. 22. Convergence curve when β is negative.
By solving this equation, the geocentric trajectory was determined. Meanwhile, Newton-Raphson iteration algorithm was applied to correct the parameters to make the whole trajectory time continuous. In the end, three example missions were established to validate the effectiveness of the design method. Declaration of competing interest The authors declare that they have no competing financial interests. Acknowledgements This work was supported by the National Natural Science Foundation of China (61833009, 61690212, 91438202, 11972130), and the National Key Research and Development Plan (grant number 2016YFB0500901). References Fig. 21. Convergence curve of the second case when β is positive.
was used to calculate the parameters of selenocentric trajectory. In order to solve the orbital problem that starts from a certain point to an excess velocity vector, two auxiliary frame system were established and a sixth-degree closed form equation was obtained.
[1] R. Yuan, C. Pingyuan, L. Enjie, Interplanetary optimum two-impulse transfer trajectories, Acta Aeronaut. Astronaut. Sin. 28 (6) (2007) 1307–1311. [2] D. Qiao, P. Cui, H. Cui, The Nereus asteroid rendezvous mission design with multi-objective hybrid optimization, in: International Symposium on Systems & Control in Aerospace & Astronautics, IEEE, 2006. [3] C.Y. Li, W.X. Jing, C.S. Gao, Adaptive backstepping-based flight control system using integral filters, Aerosp. Sci. Technol. 13 (2–3) (2009) 105–113.
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.13 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
[4] M. Yu, C. Li, Robust adaptive iterative learning control for discrete-time nonlinear systems with time-iteration-varying parameters, IEEE Trans. Syst. Man Cybern. Syst. (2017) 1–9. [5] J. Li, C. Gao, C. Li, et al., A survey on moving mass control technology, Aerosp. Sci. Technol. 82–83 (2018) 594–606. [6] T.N. Edelbaum, Optimal nonplanar escape from circular orbits, AIAA J. 9 (12) (1971) 2432–2436. [7] D. Jones, C. Ocampo, Optimal impulsive escape trajectories from a circular orbit to a hyperbolic excess velocity vector, in: American Institute of Aeronautics and Astronautics AIAA/AAS Astrodynamics Specialist Conference, Toronto, Ontario, Canada, 2010. [8] D.R. Jones, C. Ocampo, Optimization of impulsive trajectories from a circular orbit to an excess velocity vector, J. Guid. Control Dyn. 35 (1) (2012) 234–244. [9] R.J. Gerbracht, P.A. Penzo, Optimum three-impulse transfer between an elliptic orbit and a non-coplanar escape asymptote, in: AAS/AIAA Astrodynamics Specialist Conference, 1968 Sep 3-5. [10] G. Zhang, X. Zhang, X. Cao, Tangent-impulse transfer from elliptic orbit to an excess velocity vector, Chin. J. Aeronaut. 27 (3) (2014) 577–583. [11] D.H. Cho, D. Lee, H. Bang, et al., B-plane targeting method for orbit maneuver using low thrust, Int. J. Control. Autom. Syst. 15 (6120) (2017) 1–9. [12] J. Li, S. Gong, X. Wang, Analytical design methods for determining Moon-toEarth trajectories, Aerosp. Sci. Technol. 40 (2015) 138–149. [13] Yu-zhu Bai, Yu-dong Gao, Xiao-ning Xi, Fast searching design method for return transfer trajectory of lunar probe, J. Astronaut. 29 (3) (2008) 765–771. [14] Z. Lei, Design of moon return trajectory with direct atmospheric reentry, Spacecr. Eng. (2010). [15] Shengping G, Junfeng L I, Jingyang L I, Characteristics of Manned Lunar TransEarth Trajectories, Manned Spaceflight, 2012.
13
[16] S. Hong-Xin, L.I. Hai-Yang, P. Qi-Bo, et al., Point return trajectory characteristics analysis for a lunar spacecraft, J. Nat. Univ. Defense Technol. (2011). [17] Q. Luo, J. Yin, C. Han, Design of Earth—Moon free-return trajectories, J. Guid. Control Dyn. 36 (1) (2013) 263–271. [18] J. Li, S. Gong, H. Baoyin, Generation of multi-segment lunar free-return trajectories, J. Guid. Control Dyn. 36 (3) (2013) 765–775. [19] H. Wen-De, X.I. Xiao-Ning, W. Wei, et al., Characteristic analysis and design of free return orbit for lunar manned landing based on the double two-body model, J. Astronaut. (2010). [20] M. Jesick, C. Ocampo, Automated generation of symmetric lunar free-return trajectories, J. Guid. Control Dyn. 34 (1) (2011) 98–106. [21] L. Xinglong, D. Guangren, Designs of optimal free-return orbit for moon landing, in: IEEE International Conference on Robotics & Biomimetics, IEEE, 2008. [22] J. Kauffman, A successful failure: NASA’s crisis communications regarding Apollo 13, Public Relat. Rev. 27 (4) (2001) 437–448. [23] S. Belz, B. Ganzer, E. Messerschmid, et al., Hybrid life support systems with integrated fuel cells and photobioreactors for a lunar base, Aerosp. Sci. Technol. 24 (1) (2013) 169–176. [24] B. Yang, H. Yang, S. Li, Hierarchical approach for fast searching optimal launch opportunity in a wide range, Adv. Space Res. (2018). [25] Yu-zhu Bai, Xiao-ning Xi, Lei Liu, et al., Characteristics of return trajectory of lunar probe, J. Nat. Univ. Defense Technol. 30 (4) (2008) 11–16. [26] G. Potential, Orbital Mechanics for Engineering Students, Elsevier Aerospace Engineering, 2014, pp. 69–74. [27] M. Lv, M. Tan, D. Zhou, Design of two-impulse Earth–Moon transfers using differential correction approach, Aerosp. Sci. Technol. 60 (2017) 183–192, S1270963816310616.
[m5G; v1.261; Prn:7/11/2019; 15:22] P.1 (1-13)
Aerospace Science and Technology ••• (••••) ••••••
Contents lists available at ScienceDirect
Aerospace Science and Technology www.elsevier.com/locate/aescte
Deep space exploration orbit design departing from circumlunar orbit of lunar base Shi Qiu a , Xibin Cao a , Feng Wang a,∗ , Chengfei Yue a , Zhe Zhang b a b
Research Center of Satellite Technology, Harbin Institute of Technology, Harbin 150001, China Lunar Exploration and Space Engineering Center, China
a r t i c l e
i n f o
Article history: Received 15 May 2019 Received in revised form 11 October 2019 Accepted 18 October 2019 Available online xxxx Keywords: Circumlunar orbit Lunar base Fixed excess velocity vector Analytical trajectory design Genetic algorithm
a b s t r a c t A design method for a new kind of deep space mission which departs from a circumlunar orbit of the lunar base is presented in this paper. The difference between this kind of deep space mission and the traditional one is well analyzed. The method of patched conics is employed to divide the mission trajectory into four independent parts. The independent selenocentric segment is designed using a Lambert Problem based method, while the geocentric trajectory from a certain point in threedimensional space to a fixed excess velocity vector is determined analytically. Then the combination of these segments is optimized via genetic algorithm by minimizing the total velocity increasement. Finally, typically numerical simulations are demonstrated to verify the effectiveness of the proposed method. © 2019 Elsevier Masson SAS. All rights reserved.
1. Introduction Deep space exploration has significant importance in exploring the distant regions of outer space. In deep space exploration, orbit design and the orbital optimization [1,2] in various scenarios are usually the first and one of the most important steps before the spacecraft functions well under virous control strategy [3–5]. For traditional deep space missions departing from a low circular geocentric parking orbit, the orbit design problem and optimal maneuver has been studied intensively [6], and feasible impulsive transfers to some fixed hyperbolic excess velocity vector have been obtained [7,8]. The obtained impulsive solutions can be used as the initial values to determine finite-burn or continuous transfers. Among these impulsive solutions, the three-impulse time-free initial guess method reflects the behavior of the optima, resulting in better overall convergence, while the three-impulse time-fixed one and optimization method allows for quick cost assessments for a variety of escape geometries as a function of the total trip time [9,10]. However, these impulsive based methods generate orbital error in the practical application, which should be adjusted by trajectory correction maneuver using some strategies, for instance, the B-plane targeting method [11].
*
Corresponding author. E-mail address: [email protected] (F. Wang).
https://doi.org/10.1016/j.ast.2019.105505 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.
For the deep space exploration, lunar exploration and the exploration from lunar base has attracted a great deal of attention. Many countries have established different exploration missions, such as the Apollo Missions, Chang’e Missions. In these missions, the earth-to-moon transfer orbit and the return orbit design is the preliminary task. Many scholars have studied the earth-to-moon transfer orbits and the return orbits, and the design results have been verified by practical missions [12–16]. For some unexpected conditions such as a brake failure at the perigee of the moon, the free return orbits offers an optional way of coming back [17–21], and the application in Appollo-13 was called “a successful failure” [22]. Based on these remarkable results, regenerative and sustainable hybrid life support systems for a lunar base has even been designed [23]. All these improvements in lunar exploration pave the way of deep space exploration departs from lunar base. In this paper, a design method for a new kind of deep space mission which departs from a circumlunar orbit of lunar base to another planet is presented. Firstly, the difference to a traditional deep space mission starts from a low geocentric orbit is analyzed in section 2. Following that in section 3, the selenocentric trajectory is calculated on the basis of Lambert problem. The formation of initial position and velocity vectors of the geocentric trajectory are also obtained. In section 4, a new orbit design problem that departs from a certain point to a known excess velocity is presented. In order to solve this problem, two auxiliary frame system were established, thus analytic solutions are derived from a sixth-degree
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.2 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
2
equation. Next, Newton-Raphson iteration algorithm is applied to make the whole trajectory time-continuous and the genetic algorithm is used to find the best solution in section 5. Examples of deep space missions from a circumlunar orbit are demonstrated to validate the effectiveness in section 6. Finally, section 7 provides concluding remarks of this paper. 2. Problem description Launch opportunity determination is the first step of a deep space orbit design. When the ranges of launch time and flight time are settled, the best launch opportunity with the smallest total velocity increasement during the interplanetary transferring can be found using pork-chop method [24]. Once the launch opportunity is settled, the excess velocity vector denoted by v ∞ and the corresponding time instant t ∞ can be determined as well. A traditional deep space exploration mission usually starts from a low geocentric parking orbit. In these missions, when the excess velocity, the height of the acceleration point and orbital inclination are determined, the B-plane method can be applied to calculate the position and the velocity of the acceleration point. Then the geocentric hyperbolic orbit from the acceleration point to the excess velocity is determined. For the mission starting from a circumlunar orbit, the design method is quite different. The patched conics method is employed to divide the mission trajectory into four parts: the departure trajectory relative to the moon; the hyperbolic orbit relative to the earth, the cruise ellipse relative to the sun, and the hyperbolic arrival trajectory relative to the target planet. Firstly, the spacecraft is mainly influenced by the gravity of the moon, and it has to get rid of the lunar gravity in order to get into the sphere of influence of the earth, as shown in Fig. 1. The connection point of the selenocentric segment and the geocentric segment is the exit point through which the spacecraft flies out of the sphere of influence of the moon. This exit point is naturally the beginning of the geocentric hyperbolic orbit. However, this exit point can be arbitrary on the lunar sphere of influence, which leads to a geocentric orbit design problem starting from a certain point rather than a parking orbit to an excess velocity. Once this geocentric hyperbolic orbit with the excess velocity vector v ∞ is determined, the rest parts of this mission is similar to a traditional mission. Thus, the last two parts of the orbit in this mission will not be discussed in detail in this article. So, the purpose of this paper is to find an optimal orbit with a fixed excess velocity optimized over the exit point on the influence sphere of the moon as illustrated in Fig. 1. 3. Selenocentric segment design method In order to combine the selenocentric segment with the geocentric segment easily, the moon orbital plane frame coordinate O L -xyz is established [25]. The origin O L is the center of the moon, and the x- y plane is the moon’s orbital plane at the time instant t B when the spacecraft arrives at the boundary of the lunar sphere of influence. The x axis directs from the moon towards the earth while the z axis is normal to the moon’s orbital angular momentum vector. In addition, the latitude of the lunar orbit is set to zero while x-z plane is the meridian plane as is shown in Fig. 2. The point where the spacecraft flies out of the sphere of the lunar influence is denoted by B. Assume that the latitude of B is ϕ B while the longitude is λ B . Thus, the position vector of exit B can be expressed as
⎡
⎤
cos ϕ B · cos λ B L r BL = ⎣ cos ϕ B · sin λ B ⎦ · ρsoi sin ϕ B
(1)
Fig. 1. Deep space mission from a selenocentric orbit.
Fig. 2. Moon orbital plane frame coordinate. L where ρsoi is the radius of sphere of influence of the moon, and the specific value is 66200 kilometers. In this mission all six initial orbital elements are set, then the position information of the spacecraft on the parking orbit can be acquired easily at any given time with the help of high-precision gravitational model. In addition, the position of exit B can be obtained by equation (1). Once the trip time from the parking orbit to the exit is acquired, by solving the Lambert problem about T L (which will be discussed in section 5), the orbit of selenocentric segment is determined. According to the definition of the moon orbital plane frame coordinate O L -xyz, it is easy to calculate its orthogonal unit basis vectors with respect to the earth-J2000 system as follows
L i L = r E , |r EL |
k L =
r EL × v LE , |r EL × v LE |
j L = k L × i L
(2)
where, r EL and v LE are the moon’s position and the velocity relative to the earth at time instant t B when the spacecraft arrives at the exit of lunar sphere of influence, and they can be both obtained by searching high-precision ephemeris, DE421 for example. Thus, the transformation from moon orbital plane frame coordinate O L -xyz to the earth-J2000 system yields:
R LE = [i L , j L , k L ]
(3)
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.3 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
3
Fig. 4. Hyperbolic orbit description frame.
√ sin |β| =
Fig. 3. Hyperbolic orbit description frame.
Therefore, the end position of the selenocentric segment which is the initial position of the geocentric segment relative to the earth can be expressed as follows
r BE = R LE · r BL + r LE
(4)
4. Geocentric segment design method Till now, the initial position r BE and hyperbolic excess velocity v ∞ with respect to the earth-J2000 frame are obtained. Then the geocentric orbit design problem from a certain position r BE to an excess velocity vector v ∞ will be discussed in this section. As shown in Fig. 3, the earth lies at one of the focus of hyperbolic orbit. In order to describe the orbit easily, we established a new reference frame H-xyz. The x- y plane of the H-xyz frame is the plane of the orbit, and the x axis is directed along the apse line, and the center H is |a| + r p away from the earth. The z axis is normal to the plane of the orbit in the direction of the orbital angular momentum vector. In this new frame, the orbit can be described in hyperbolic equation as
x2
|a|2
y2
−
|b|2
=1
(5)
where a is the semimajor axis while b is the semiminor axis that can be found in terms of a and the eccentricity e of the hyperbolic orbit as
b = a e2 − 1
(6)
According to [26], the specific energy of a hyperbolic orbit is
ε=
μE
(7)
2|a|
with μ E the gravitational parameter of the earth. Substituting ε , v ∞ into to living force equation yields
|v ∞ |2 2 where
|a| =
−
μE μE = E 2|a| ρsoi
(8)
E ρsoi is the radius of sphere of influence of the earth. Thus,
μE |v ∞ |2 − 2ρμE E
(9)
soi
The angle β in Fig. 3 is the angle between asymptote and apse, which satisfies
e2 − 1
cos |β| =
,
e
1
(10)
e
To get the transformation from the earth-J2000 frame to H-xyz system, an auxiliary frame system M-xyz is established. The x axis of M-xyz is directed along the asymptote of the hyperbola and the z axis is normal to the plane of the orbit in the direction of the orbital angular momentum vector, as is shown in Fig. 4. Considering that the excess velocity is parallel to the asymptote when they approach infinity, the orthogonal unit basis vectors for this auxiliary frame can be expressed as follows
i M = v ∞ , |v ∞ |
k M =
r0 × v ∞ , |r0 × v ∞ |
j M = k M × i M
(11)
Then the transformation from the earth-J2000 frame to the Mxyz system can be obtained as:
T T T RM J = i M , j M , kM
(12)
Then the initial position vector r BE in the M-xyz frame system is
r M = R MJ · r BE
(13)
By now, it is easy to obtain the transformation from M-xyz to H-xyz by rotating a angle β about the z axis. The sign of angle β is positive when the rotation is clockwise. Otherwise, the angle β is negative. Thus, transformation from M-xyz to H-xyz is
⎡
H RM
cos β = ⎣ − sin β 0
⎤
sin β cos β 0
0 0⎦ 1
⎧⎡ e 2 −1 1 ⎪ ⎪ e e ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ e 2 −1 1 ⎪ ⎣− e ⎪ ⎪ e ⎪ ⎪ ⎨ 0 0 = ⎡ 2 ⎪ 1 ⎪ ⎪ − e e −1 ⎪ e ⎪ ⎢ ⎪ ⎪ ⎢ 2 ⎪ 1 ⎪ ⎣ e −1 ⎪ e e ⎪ ⎩
0
0
⎤ 0
⎥
⎥ 0⎦ β >0 1 0
⎤
(14)
⎥
⎥ 0⎦ β <0 1
Therefore, the transformation from the earth-J2000 to H-xyz can be calculated by (12) and (14): H R EH = R M · RM E
(15)
Furthermore, the origin of H-xyz system is |a| + r p behind the origin of the earth-J2000 system. The initial position r BE relevant the H-xyz system can be expressed as follows
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.4 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
4
⎡ H r H = R M · RM r BE + ⎣ E ·
|a| + r p 0 0
⎤
⎡
|a|e
⎤
H ⎦ = RM · r M + ⎣ 0 ⎦
(16)
0
where r p is the perigee radius of the hyperbolic orbit that can be found in terms of a and e as
r p = |a|(e − 1) Substituting
(17) H RM
into to equation (16) yields
⎧⎡ ⎤ e 2 −1 1 M ⎪ ⎪ rx + e r M + |a|e y ⎪ e ⎪ ⎥ ⎢ ⎪ ⎪ ⎥ β >0 ⎢ ⎪ e 2 −1 M 1 M ⎪ ⎦ ⎣ − r + r ⎪ x e e y ⎪ ⎪ ⎪ ⎨ 0 r H = ⎡ ⎤ 2 ⎪ e −1 M 1 M ⎪ ⎪ r − r + | a | e y ⎪ e ⎪ ⎢ e x ⎥ ⎪ ⎪ ⎢ ⎥ β <0 ⎪ e 2 −1 M ⎪ ⎣ ⎦ r x + 1e r M ⎪ y e ⎪ ⎪ ⎩
(18)
0
Fig. 5. Example of multiple solutions and restricted area.
where r xM and r M r M . The hypery are the first two components of bolic orbit lies in the x- y plane of the H-xyz system, so the third component of r M is zero. Taking the case β > 0 as example, substituting r H into to the hyperbolic equation, we can obtain
( 1e r xM +
e 2 −1 M ry e a2
Denoting x =
+ |a|e )2
−
(−
√
e 2 − 1 yields
2 r xM + r M y x + |a| x − 1
2
e 2 −1 M 2 r x + 1e r M y ) e 2 2 a (e − 1)
x2 − −r xM x + r M y
2
=1
(19)
= |a|x2 x2 − 1 (20)
The above equation can further be written in the following formation. 6
5
4
3
2
c6 x + c5 x + c4 x + c3 x + c2 x + c1 x + c0 = 0
(21)
Therefore, a sixth-degree closed form equation with only one variable x is obtained. All the coefficients of this equation are expressed as follows:
⎧ ⎪ c = a2 ⎪ ⎪ 6 ⎪ ⎪ ⎪ c 5 = 2|a|r M ⎪ y ⎪ ⎪ ⎪ M 2 ⎪ ⎪ ⎪ c = ry + 2|a|r xM + a2 ⎪ ⎨ 4 M c 3 = 2r xM r M y + 2|a|r y ⎪ ⎪ ⎪ c = 2|a|r M ⎪ ⎪ 2 x ⎪ ⎪ ⎪ M M ⎪ ⎪ c 1 = 2r x r y ⎪ ⎪ ⎪ ⎪ ⎩ c = −r M 2 0
r x2
−
a2
r 2y a2 (e 2 − 1)
H r Bx
a2
x−
y
(23)
(24)
r BHy a2 (e 2 − 1)
y=1
(25)
Based on the energy conservation principle, the magnitude of velocity at exit B can be obtained from the following equation.
1 2
v 2B −
μE |r B |
=
1 2
v 2∞ −
μE |r∞ |
(26)
Then the velocity | v B | can be obtained:
|v B | =
(22)
=1
According to tangent formula of hyperbolic, the direction of velocity at exit B in geocentric segment must be parallel to the following line:
Similarly, when β is negative, we can also acquire a closed form equation of sixth degree. The coefficients are
⎧ ⎪ c 6 = a2 ⎪ ⎪ ⎪ ⎪ ⎪ c 5 = −2|a|r M ⎪ y ⎪ ⎪ ⎪ M 2 ⎪ ⎪ ⎪ c = ry + 2|a|r xM + a2 ⎪ ⎨ 4 M c 3 = −2r xM r M y − 2|a|r y ⎪ ⎪ ⎪ ⎪ c 2 = 2|a|r xM ⎪ ⎪ ⎪ ⎪ M M ⎪ ⎪ ⎪ c 1 = −2r x r y ⎪ ⎪ ⎪ ⎩ c = −r M 2 0 y
By solving the equation, the eccentricity of the hyperbola can be obtained, so that the hyperbolic orbit equation in the H-xyz system yields
v 2∞ +
2μ E
|r B |
−
2μ E
|r∞ |
(27)
It is worth noticing that the direction of the orbital angular momentum is relevant to flight direction of the orbit. Therefore, there may be more than one orbits, one clockwise hyperbolic orbit and an anticlockwise one, departing from a certain position to a given excess velocity vector satisfying the boundary conditions. What’s more, for each certain excess velocity v ∞ , there is a limited area where only one hyperbolic orbit can be found satisfying given conditions. Take Fig. 5 as an example, the left boundary of the shadow region is parallel to the excess velocity vector while the right one is symmetric with the left one about x axis. Therefore, the shadow region will be feasible region of the anticlockwise hyperbolic obits with the earth as their focus and the exclusion region of the clockwise orbits. However, solving the equation (21) may result in some infeasible solutions because the excess velocity vector only represents the direction not the position when the spacecraft flies out of the influence sphere of the earth. For example, the hyperbolic orbit 3 in Fig. 5 is an infeasible solution passing through the point P and with a same velocity magnitude of v ∞ . In order to reject these infeasible solutions, we can calculate the angle γ between the excess
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.5 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
5
In order to make the entire trajectory time continuous, the time t B has to be chosen consistently such that: t∞ = t∞
(31)
Then t B can be chosen via Newton-Raphson iteration according to [27]. Up to now, the geocentric trajectory along with the heliocentric one has been determined. According to the analysis of section 1, once the trip time in the lunar sphere of influence is settled, the selenocentric trajectory can be acquired by solving the Lambert problem. Thus, the departure velocity and the velocity of the arrival at the exit are both obtained. Therefore, the total velocity increasement yields
− v = |v d − v 0 | + v + B − vB
Fig. 6. Example of virtual branch of hyperbolic orbit.
velocity vector v ∞ and the position vector r EA , in which E represents the earth while A represents the exit point of the earth’s influence sphere.
γ = arc cos
v ∞ · r EA
(28)
If γ is an acute angle, this solution is a feasible one. Otherwise when γ is an obtuse angle, the solution will be infeasible. Considering the symmetricity of the hyperbola, the real trajectory is always accompanied with its reflection. Mathematically, if we replace x and y with −x and − y, we will obtain the same equation:
(−x)2 (− y )2 x2 y2 − = − =1 |a|2 |b|2 |a|2 |b|2
(29)
Thus, the solution to the equation (21) may also result in a virtual trajectory as is shown in Fig. 6, which forms the infeasible hyperbolic orbit 3 in Fig. 5. 5. Optimization of the orbit When the pork-chop method is applied to finding the best opportunities to travel from the earth to another planet, the time t ∞ of the excess velocity v ∞ is determined. Meanwhile, the solution to the geocentric segment is only depend on the initial position vector r BE and the excess velocity v ∞ , which means the trip time of the geocentric segment denoted by T E is fixed when r BE and v ∞ are chosen. For the initial position vector r BE , it is determined by the latitude and longitude of the exit of lunar sphere of influence as well as the moon’s real-time position vector r LE in the sphere of influence of the earth. Hence, if the time t B when the spacecraft arrives at the exit of lunar sphere of influence is determined, the real time that the spacecraft flies out of the sphere of influence of the earth is settled as well. t∞ = tB + T E
where, v d represents the velocity of departure while v 0 is the velocity of the parking orbit. v + B is the initial velocity of the geocentric trajectory while v − B is the ending velocity of the selenocentric trajectory. Since the total velocity increasement is not an analytical solution to the trip time of the selenocentric segment T L , an intelligent optimization algorithm can be employed to minimize the total velocity increasement with respect to T L . Then the optimization index, the total velocity increasement, can be expressed by a function of T L as:
− J = v = | v d − v 0 | + v + B − v B = f (T L )
|v ∞ | · |r EA |
(32)
(30)
where, t ∞ is the real time that the spacecraft flies out of the sphere of influence of the earth.
(33)
To obtain the optimal solution of equation (33) the genetic algorithm, owing to its significant advantage of the capability of initial-condition-independent global optimization and the strong robustness, has been adopted. The procedure of solving this optimization problem is listed as follows and Fig. 7. Step 1. Using pork-chop method to find the best opportunity departing from the earth to another planet, i.e., to find the excess velocity vector v ∞ and its time t ∞ ; Step 2. By choosing the latitude λ B and longitude ϕ B of the exit of lunar sphere of influence, as well as the time t B , to get the initial position relative to the earth r BE ; Step 3. Using the results of v ∞ and r BE and the method introduced in section 3 to get geocentric segment; Step 4. Calculate the trip time from the initial position to the excess velocity vector v ∞ ; and verify whether it Step 5. Calculate the real arrived time t ∞ satisfies the time error condition or not; Step 6. Repeat Step 2 to Step 5 until the real arrived time t ∞ satisfies the time error condition; Step 7. Calculate the optimal total velocity increasement. 6. Simulation cases In order to validate the effectiveness of the design method, three simulation cases with different target planets are presented in this section. For the first two examples, Mars exploration missions are considered while the third one considers the Venus exploration mission. Thus, the effectiveness of proposed method with different launching time and different target planets are both verified through the listed scenarios. Example 1 Assume a deep space exploration mission to the mars which will be carried out from a circumlunar orbit during 2020. The initial orbital elements relevant to the moon − J 2000 frame of a circumlunar spacecraft are shown in Table 1.
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.6 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
6
Fig. 7. Flow chart of the mission orbit design. Table 1 Initial orbital elements of a circumlunar orbit. Time (UTCG)
Semimajor axis (km)
Eccentricity
Inclination (deg)
RAAN (deg)
Arg of perigee (deg)
True anomaly (deg)
1 Jul 2020 00:00:00.000
8943.996897
0.693000
25.400
275.800
182.400
341.200
Fig. 9. Three-dimensional figure of total velocity increasement in 2020. Fig. 8. Pork-chop figure of total velocity increasement in 2020.
According to the analysis in this article, the excess velocity vector v ∞ and its time t ∞ have to be calculated at first. A pork-chop figure about total velocity increasement is shown in Fig. 8 (Fig. 9). It is easy to find the best opportunity using an optimal method in this pork-chop figure, marked with a star symbol “∗”. The
departure time and the excess velocity of the best opportunity yield
t ∞ = 2459054.7534 (Julian day) v ∞ = [3.2298; 0.7070; 1.6406] (km/s)
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.7 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
7
Table 2 Coefficients of equation (21) when β is positive. Coefficients
c 0 (e+11)
c 1 (e+10)
c 2 (e+9)
c 3 (e+9)
c 4 (e+11)
c 5 (e+10)
c 6 (e+8)
Values
-1.1620
1.1752
-1.0769
-9.5428
1.1610
-2.1295
9.7565
Table 3 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
11.8796
9.8612
1.0521
0.0+1.0i
0.0-1.0i
-0.9663
Table 4 Coefficients of equation (21) when β is negative. Coefficients
c 0 (e+11)
c 1 (e+10)
c 2 (e+9)
c 3 (e+9)
c 4 (e+11)
c 5 (e+10)
c 6 (e+8)
Values
-1.1620
-1.1752
-1.0769
9.5428
1.1610
2.1295
9.7565
Table 5 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
-11.8796
-9.8612
-1.0521
0.0+1.0i
0.0-1.0i
0.9663
Table 6 Orbital elements of feasible solutions of geocentric segment. Conditions
Julian day of exit B
Eccentricity of geocentric segment
r B in earth-J2000 (km)
v B in earth-J2000 (km/s)
β >0
2459052.0256 2459052.0369
11.7708 9.7677
[−169592.6; 257181.3; 135567.6] [−170361.7; 256529.5; 135359.7]
[3.2529; 1.0112; 1.8777] [3.229; 1.057; 1.8995]
β <0
2459051.1055
1.3070
[−101775.2; 303652.0; 149094.7]
[0.8726; −3.3514; −1.7391]
To a common deep space mission, it usually takes about 3 days to fly out of the sphere of influence of the earth. Thus, the initial guessed value of t B can be assumed to be 3 days earlier than t ∞ . Meanwhile, in this example case, the latitude and longitude of the exit of lunar sphere of influence are both set to zero. When β is positive, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 2. It is easy to obtain the solutions to the six-degree equation which are listed in Table 3. It is apparently that x4 , x5 and x6 are infeasible solutions because the eccentricity have to be a real value and must be greater than 1. Therefore, potential feasible solutions must be among x1 , x2 and x3 . As for the solution x1 = 11.8796, initial trip time is 2.6581 days, not equal to 3 days. Using Newton-Raphson iteration algorithm, we can obtain convergent value, x1 = 11.7282, under which condition, the trip time of geocentric segment is 2.7278 days. Similarly, for the solution x2 = 9.8612, the convergent value is x2 = 9.7164 and the geocentric trip time is 2.7165 days. However, for the solution x3 = 1.0521, the excess velocity vector v ∞ relative to H-xyz frame is [2.5428; −2.6753; 0.0000] while the velocity when the spacecraft flies out of the influence sphere of the earth is [−2.5428; −2.6753; 0.0000]. That means this solution lies on the other branch of the hyperbolic orbit which does not exit indeed. When β is negative, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 4. It is easy to obtain the solutions to the six-degree equation which are listed in Table 5. It is obviously that only x6 may be a feasible solution. The convergent value is x1 = 0.8415 and the geocentric trip time is 3.6479 days.
Fig. 10. Convergence curve of the first case when β is positive.
The orbital elements of feasible solutions are listed in Table 6. Thus, the geocentric trajectory has been obtained. The next step is to find a best selenocentric solution to minimize the total velocity increasement the least according to the optimal method in section 5. Genetic algorithm is applied in this case. The convergence curves under different conditions are shown in Figs. 10 to 12. It is obvious that genetic algorithm is suitable in this optimal problem for its fast convergence. The best solutions are illustrated in Table 7.
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.8 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
8
Table 7 Orbital elements of selenocentric segment. Conditions
Selenocentric trip time
Acceleration position in moon-J2000 (km)
Acceleration velocity in moon-J2000 (km/s)
β >0
1.9315 days 2.8296 days
[−4636.9; −8371.5; −2592.2] [−4748.4; −7881.5; −2621.3]
[0.2788; −0.8237; −0.3637] [0.1750; −0.8481; −0.3698]
β <0
0.5507 days
[5639.9; −12133.7; 2082.1]
[−0.7089; 1.5064; −0.1930]
Fig. 13. Pork-chop figure of total velocity increasement in 2022.
Fig. 11. Convergence curve of the second case when β is positive.
Fig. 12. Convergence curve when β is negative.
Example 2 Assume another deep space exploration mission to the mars which will be carried out from a circumlunar orbit during 2021. The initial orbital elements relevant to the moon − J 2000 frame of a circumlunar spacecraft are shown in Table 8.
Fig. 14. Three-dimensional figure of total velocity increasement in 2022.
The pork-chop figure about total energy consumption in this case is shown in Fig. 13 (Fig. 14).
Table 8 Initial orbital elements of a circumlunar orbit. Time (UTCG)
Semimajor axis (km)
Eccentricity
Inclination (deg)
RAAN (deg)
Arg of perigee (deg)
True anomaly (deg)
1 Sep 2022 00:00:00.000
8943.996897
0.693000
25.400
275.800
182.400
341.200
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.9 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
9
Table 9 Coefficients of equation (21) when β is positive. Coefficients
c 0 (e+10)
c 1 (e+11)
c 2 (e+10)
c 3 (e+11)
c 4 (e+10)
c 5 (e+9)
c 6 (e+8)
Values
-3.8284
-1.4887
1.4648
-1.5640
5.3303
-7.5338
3.7064
Table 10 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
10.3998
5.0816+4.0197i
5.0816-4.0197i
0.0+1.0i
0.0-1.0i
-0.2366
Table 11 Coefficients of equation (21) when β is negative. Coefficients
c 0 (e+10)
c 1 (e+11)
c 2 (e+10)
c 3 (e+11)
c 4 (e+10)
c 5 (e+9)
c 6 (e+8)
Values
-3.8284
1.4887
1.4648
1.5640
5.3303
7.5338
3.7064
Table 12 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
-10.3998
5.0816+4.0197i
5.0816-4.0197i
0.0+1.0i
0.0-1.0i
0.2366
Table 13 Orbital elements of feasible solutions of geocentric segment. Conditions
Julian day of exit B
Eccentricity of geocentric segment
r B in earth-J2000 (km)
v B in earth-J2000 (km/s)
β >0 β <0
2459840.6552 2459838.8930
12.2782 1.0274
[83290.5; 359056.7; 182649.6] [225503.9; 329137.8; 157305.9]
[2.4096; 2.2604; 3.4244] [−2.5021; −3.6221; −1.7841]
It is easy to find the best opportunity using an optimal method in this pork-chop figure, marked with a star symbol “∗”. The time and the excess velocity of the best opportunity yield
t ∞ = 2459530.0293 v ∞ = [0.3589; −4.1368; 0.7600] Also, in this example, the initial guessed value of t B is assumed to be 3 days earlier than t ∞ and the latitude and longitude of the exit of lunar sphere of influence are both set to zero. When β is positive, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 9. It is easy to obtain the solutions to the six-degree equation which are listed in Table 10. Obviously, only x1 can be a potential feasible solution. The convergent value is x1 = 12.2372 and the geocentric trip time is 1.3679 days. When β is negative, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 11. It is easy to obtain the solutions to the six-degree equation which are listed in Table 12. It is obviously that only x6 may be a feasible solution. The convergent value is x6 = 0.2357 and the geocentric trip time is 3.1301 days. The orbital elements of feasible solutions are listed in Table 13. Thus, the geocentric trajectory has been obtained. The next step is to find a best selenocentric solution to make the total velocity increasement the least according to the optimal method in section 5. Genetic algorithm is applied in this case. The convergence curves under different conditions are shown in Fig. 15 and Fig. 16.
Fig. 15. Convergence curve when β is positive.
It is obvious that genetic algorithm is suitable in this optimal problem for its fast convergence speed. The best solutions are displayed in Table 14. Nevertheless, it is worth noticing that, in the case β is negative, the spacecraft will fly by the earth rather than directly fly out of the sphere of influence of the earth when β is positive. Meanwhile, the perigee radius of the hyperbolic orbit in the negative case is smaller than the radius of the earth. That means, the spacecraft will intersect into the earth, as shown in Fig. 17. Therefore, the trajectory in the negative case cannot be applied into engineering.
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.10 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
10
Example 3 Assume another deep space exploration mission to the Venus which will be carried out from a circumlunar orbit during 2023. The initial orbital elements relevant to the moon − J 2000 frame of a circumlunar spacecraft are shown in Table 15. The pork-chop figure about total energy consumption in this case is shown in Fig. 18 (Fig. 19). It is easy to find the best opportunity using an optimal method in this pork-chop figure, marked with a star symbol “∗”. The time and the excess velocity of the best opportunity yield
t ∞ = 2459530.0291 v ∞ = [−1.8218; 2.8860; −0.6494]
Fig. 16. Convergence curve when β is negative.
Also, in this example, the initial guessed value of t B is assumed to be 3 days earlier than t ∞ and the latitude and longitude of the exit of lunar sphere of influence are both set to zero. When β is positive, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 16. It is easy to obtain the solutions to the six-degree equation which are listed in Table 17.
Table 14 Orbital elements of selenocentric segment. Conditions
Selenocentric trip time
Acceleration position in moon-J2000 (km)
Acceleration velocity in moon-J2000 (km/s)
β >0 β <0
3.1925 days 0.5507 days
[6215.8; −3382.3; 2774.1] [6186.2; −3282.0; 2764.9]
[−0.5033; 0.9420; 0.0928] [−1.2829; 1.1034; −0.4854]
Fig. 17. Trajectory in earth-J2000 frame when β is negative.
Table 15 Initial orbital elements of a circumlunar orbit. Time (UTCG)
Semimajor axis (km)
Eccentricity
Inclination (deg)
RAAN (deg)
Arg of perigee (deg)
True anomaly (deg)
25 May 2023 00:00:00.000
2037.400338
0
45.398
32.335
190.077
43.793
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.11 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
11
It is apparently that x3 , x5 and x6 are infeasible solutions because the eccentricity have to be a real value and must be greater than 1. Therefore, potential feasible solutions must be among x1 , x2 and x4 . As for the solution x1 = 10.1589, initial trip time is 3.2148 days, not equal to 3 days. Using Newton-Raphson iteration algorithm, we can obtain convergent value, x1 = 10.2052, under which condition, the trip time of geocentric segment is 3.1749 days. Similarly, for the solution x4 = 0.7054, the convergent value is x4 = 0.8594 and the geocentric trip time is 3.8180 days. However, for the solution x2 = 8.1155, the excess velocity vector v ∞ relative to H-xyz frame is [0.4249; −3.4481; 0.0000] while the velocity when the spacecraft flies out of the influence sphere of the earth is [0.4249; 3.4481; 0.0000]. That means this solution lies on the other branch of the hyperbolic orbit which does not exit indeed. When β is negative, the initial coefficients of the six-degree equation when trip time of geocentric segment is set to 3 days are listed in Table 18. It is easy to obtain the solutions to the six-degree equation which are listed in Table 19. It is obviously that only x3 may be a feasible solution. The convergent value is x3 = 1.1401 and the geocentric trip time is 3.8257 days. The orbital elements of feasible solutions are listed in Table 20. Thus, the geocentric trajectory has been obtained. The next step is to find a best selenocentric solution to make the total velocity increasement the least according to the optimal method in section 5. Genetic algorithm is applied in this case. The convergence curves under different conditions are shown in Figs. 20 to 22. It is obvious that genetic algorithm is suitable in this optimal problem for its fast convergence speed. The best solutions are displayed in Table 21.
Fig. 18. Pork-chop figure of total velocity increasement in 2023.
7. Conclusion In this paper, we introduced a new type of deep space mission departing from a circumlunar orbit. Since the initial central body is the moon, not the earth, the selenocentric trajectory is a new part of the whole mission compared with the traditional one. The mission trajectory is divided into four independent parts using patched conics method. The method based on Lambert problem
Fig. 19. Three-dimensional figure of total velocity increasement in 2023. Table 16 Coefficients of equation (21) when β is positive. Coefficients
c 0 (e+10)
c 1 (e+10)
c 2 (e+9)
c 3 (e+10)
c 4 (e+10)
c 5 (e+10)
c 6 (e+9)
Values
-9.8415
8.7773
-9.9506
6.5459
8.9729
-2.2314
1.2648
Table 17 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
10.1589
8.1155
-1.3379
0.7054
+1.0i
-1.0i
Table 18 Coefficients of equation (21) when β is negative. Coefficients
c 0 (e+10)
c 1 (e+10)
c 2 (e+9)
c 3 (e+10)
c 4 (e+10)
c 5 (e+10)
c 6 (e+9)
Values
-9.8415
-8.7773
-9.9506
-6.5459
8.9729
2.2314
1.2648
Table 19 Solutions to equation (21) when β is positive. Solutions
x1
x2
x3
x4
x5
x6
Values
-10.1589
-8.1155
1.3379
-0.7054
+1.0i
-1.0i
JID:AESCTE
AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.12 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
12
Table 20 Orbital elements of feasible solutions of geocentric segment. Conditions
Julian day of exit B
Eccentricity of geocentric segment
r B in earth-J2000 (km)
v B in earth-J2000 (km/s)
β >0
2460098.3879 2460097.7448
10.2541 1.3185
[−147868.0; −277401.4; −134052.0] [−194745.1; −248150.1; −116028.6]
[−2.2208; 2.7903; −0.9073] [2.3835; 2.4533; 1.3712]
β <0
2460097.7371
1.5165
[−195271.9; −247764.1; −115795.8]
[2.4056; 2.4296; 1.3749]
Table 21 Orbital elements of selenocentric segment. Conditions
Selenocentric trip time
Acceleration position in moon-J2000 (km)
Acceleration velocity in moon-J2000 (km/s)
β >0
1.6314 days 0.5991 days
[−615.6; −1627.2; −1060.2] [262.1; −1461.9; −1394.6]
[1.0030; −0.2403; −0.1578] [1.5493; 0.2189; 0.0565]
β <0
0.5913 days
[274.4; −1457.2; −1397.2]
[1.5619; 0.2283; 0.0609]
Fig. 20. Convergence curve of the first case when β is positive.
Fig. 22. Convergence curve when β is negative.
By solving this equation, the geocentric trajectory was determined. Meanwhile, Newton-Raphson iteration algorithm was applied to correct the parameters to make the whole trajectory time continuous. In the end, three example missions were established to validate the effectiveness of the design method. Declaration of competing interest The authors declare that they have no competing financial interests. Acknowledgements This work was supported by the National Natural Science Foundation of China (61833009, 61690212, 91438202, 11972130), and the National Key Research and Development Plan (grant number 2016YFB0500901). References Fig. 21. Convergence curve of the second case when β is positive.
was used to calculate the parameters of selenocentric trajectory. In order to solve the orbital problem that starts from a certain point to an excess velocity vector, two auxiliary frame system were established and a sixth-degree closed form equation was obtained.
[1] R. Yuan, C. Pingyuan, L. Enjie, Interplanetary optimum two-impulse transfer trajectories, Acta Aeronaut. Astronaut. Sin. 28 (6) (2007) 1307–1311. [2] D. Qiao, P. Cui, H. Cui, The Nereus asteroid rendezvous mission design with multi-objective hybrid optimization, in: International Symposium on Systems & Control in Aerospace & Astronautics, IEEE, 2006. [3] C.Y. Li, W.X. Jing, C.S. Gao, Adaptive backstepping-based flight control system using integral filters, Aerosp. Sci. Technol. 13 (2–3) (2009) 105–113.
JID:AESCTE AID:105505 /FLA
[m5G; v1.261; Prn:7/11/2019; 15:22] P.13 (1-13)
S. Qiu et al. / Aerospace Science and Technology ••• (••••) ••••••
[4] M. Yu, C. Li, Robust adaptive iterative learning control for discrete-time nonlinear systems with time-iteration-varying parameters, IEEE Trans. Syst. Man Cybern. Syst. (2017) 1–9. [5] J. Li, C. Gao, C. Li, et al., A survey on moving mass control technology, Aerosp. Sci. Technol. 82–83 (2018) 594–606. [6] T.N. Edelbaum, Optimal nonplanar escape from circular orbits, AIAA J. 9 (12) (1971) 2432–2436. [7] D. Jones, C. Ocampo, Optimal impulsive escape trajectories from a circular orbit to a hyperbolic excess velocity vector, in: American Institute of Aeronautics and Astronautics AIAA/AAS Astrodynamics Specialist Conference, Toronto, Ontario, Canada, 2010. [8] D.R. Jones, C. Ocampo, Optimization of impulsive trajectories from a circular orbit to an excess velocity vector, J. Guid. Control Dyn. 35 (1) (2012) 234–244. [9] R.J. Gerbracht, P.A. Penzo, Optimum three-impulse transfer between an elliptic orbit and a non-coplanar escape asymptote, in: AAS/AIAA Astrodynamics Specialist Conference, 1968 Sep 3-5. [10] G. Zhang, X. Zhang, X. Cao, Tangent-impulse transfer from elliptic orbit to an excess velocity vector, Chin. J. Aeronaut. 27 (3) (2014) 577–583. [11] D.H. Cho, D. Lee, H. Bang, et al., B-plane targeting method for orbit maneuver using low thrust, Int. J. Control. Autom. Syst. 15 (6120) (2017) 1–9. [12] J. Li, S. Gong, X. Wang, Analytical design methods for determining Moon-toEarth trajectories, Aerosp. Sci. Technol. 40 (2015) 138–149. [13] Yu-zhu Bai, Yu-dong Gao, Xiao-ning Xi, Fast searching design method for return transfer trajectory of lunar probe, J. Astronaut. 29 (3) (2008) 765–771. [14] Z. Lei, Design of moon return trajectory with direct atmospheric reentry, Spacecr. Eng. (2010). [15] Shengping G, Junfeng L I, Jingyang L I, Characteristics of Manned Lunar TransEarth Trajectories, Manned Spaceflight, 2012.
13
[16] S. Hong-Xin, L.I. Hai-Yang, P. Qi-Bo, et al., Point return trajectory characteristics analysis for a lunar spacecraft, J. Nat. Univ. Defense Technol. (2011). [17] Q. Luo, J. Yin, C. Han, Design of Earth—Moon free-return trajectories, J. Guid. Control Dyn. 36 (1) (2013) 263–271. [18] J. Li, S. Gong, H. Baoyin, Generation of multi-segment lunar free-return trajectories, J. Guid. Control Dyn. 36 (3) (2013) 765–775. [19] H. Wen-De, X.I. Xiao-Ning, W. Wei, et al., Characteristic analysis and design of free return orbit for lunar manned landing based on the double two-body model, J. Astronaut. (2010). [20] M. Jesick, C. Ocampo, Automated generation of symmetric lunar free-return trajectories, J. Guid. Control Dyn. 34 (1) (2011) 98–106. [21] L. Xinglong, D. Guangren, Designs of optimal free-return orbit for moon landing, in: IEEE International Conference on Robotics & Biomimetics, IEEE, 2008. [22] J. Kauffman, A successful failure: NASA’s crisis communications regarding Apollo 13, Public Relat. Rev. 27 (4) (2001) 437–448. [23] S. Belz, B. Ganzer, E. Messerschmid, et al., Hybrid life support systems with integrated fuel cells and photobioreactors for a lunar base, Aerosp. Sci. Technol. 24 (1) (2013) 169–176. [24] B. Yang, H. Yang, S. Li, Hierarchical approach for fast searching optimal launch opportunity in a wide range, Adv. Space Res. (2018). [25] Yu-zhu Bai, Xiao-ning Xi, Lei Liu, et al., Characteristics of return trajectory of lunar probe, J. Nat. Univ. Defense Technol. 30 (4) (2008) 11–16. [26] G. Potential, Orbital Mechanics for Engineering Students, Elsevier Aerospace Engineering, 2014, pp. 69–74. [27] M. Lv, M. Tan, D. Zhou, Design of two-impulse Earth–Moon transfers using differential correction approach, Aerosp. Sci. Technol. 60 (2017) 183–192, S1270963816310616.