A UKF based design method for lunar free-return trajectories

Accepted Manuscript A UKF based design method for lunar free-return trajectories

Hongli Zhang, Chao Han, Qinqin Luo, Jian Li

PII: DOI: Reference:

S1270-9638(15)00351-X http://dx.doi.org/10.1016/j.ast.2015.11.007 AESCTE 3475

To appear in:

Aerospace Science and Technology

Received date: Revised date: Accepted date:

4 August 2015 6 November 2015 9 November 2015

Please cite this article in press as: H. Zhang et al., A UKF based design method for lunar free-return trajectories, Aerosp. Sci. Technol. (2015), http://dx.doi.org/10.1016/j.ast.2015.11.007

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A UKF Based Design Method for Lunar Free-Return Trajectories Hongli Zhang* and Chao Han† Beihang University, Beijing, People’s Republic of China, 100191

Qinqin Luo‡ Beijing Aerospace Technology Institute, Beijing, People’s Republic of China, 100191

Jian Li§ Shanghai Institute of Satellite Engineering, Shanghai, People’s Republic of China, 200240

Abstract A method based on unscented Kalman filter parameter estimation is proposed for the design of lunar free-return trajectories. The method avoids calculating the Jacobian matrix and obtains large convergence ability compared with the common differential-correction method. Given that the initial estimate of the free-return trajectory is only required to be generated under the two-body Earth-spacecraft model, the difficulty of guessing a good initial estimate is greatly reduced. Through converting the original problem to a parameter estimation representation, the method is used to find the converged final solution that satisfies the constraints at injection, perilune, and Earth-entry interface under a high-fidelity gravitational model. In addition to its deceptively simplicity, the method proves to be quite effective through numerical examples in finding the solution to lunar free-return trajectories within a few iterations with great numerical accuracy.

Nomenclature Am

=

CR cl h i I Jd K Ls R P r, v r *, v * rEM

= = = = = = = = = = = = =

satellite area/mass ratio solar radiation coefficient speed of light, km/s orbital altitude, km orbital inclination, deg identity matrix Julian date Kalman gain matrix solar luminosity, kgm2/s3 radius, km state covariance matrix position and velocity in the J2000 ECI frame, km, km/s position and velocity with respect to the moon, km, km/s moon’s position in the J2000 ECI frame, km

Re Rr t

= = =

measurement noise matrix system noise matrix time, s

*

Ph.D. Candidate, School of Astronautics; [email protected]. Professor, School of Astronautics; [email protected]. ‡ Ph.D.; [email protected]. § Ph.D.; [email protected]. †

1

tEM

=

transfer time of perilune, s

TM

=

epoch of perilune, Jd

TTLI = vx , v y , vz =

epoch of Earth departure, Jd velocity components, km/s

w = W ( m) ,W (c ) = x, y , z = x = α, β = γ =

parameters weight factors for the state (m) and covariance (c) position components, km state vector scaling parameters flight-path angle, deg

μ E , μM =

gravitational parameters of Earth and moon, km3/s2 right ascension of ascending node, deg argument of perigee, deg

Ω

= = Subsripts = E M/m = o = p = re = S = = TLI = 0 = 1 2 = Superscripts = t T = ' =

ω

Earth conditions moon conditions Earth parking orbit conditions point of closest approach reentry conditions Sun conditions quantity at translunar injection initial estimate conditions relative to the primary mass conditions relative to the secondary mass target condition transpose transformed quantity

INTRODUCTION Since lunar exploration activities began in the last century, lunar free-return trajectories have drawn much attention. Lunar free-return trajectories and their variations prove useful for future human lunar missions, given that they increase the chances of crew survival in case the mission has to be aborted. The first three Apollo missions to the moon flew circumlunar free-return trajectories [1, 2]. The design of lunar free-return trajectories is essentially a two-point boundary value problem (TPBVP), which is generally divided into two steps. At first, an initial estimate is obtained by some analytical techniques. Based on the initial estimate, the final solution corresponding to the precise transfer trajectory that satisfies the constraints on targeting parameters is found using certain numerical search algorithms. For each single-flyby free-return trajectories, the time of flight from Earth departure to Earth arrival is on the order of days, whereas the time spent in the vicinity of the moon is on the order of hours. Therefore, the problem is highly sensitive to the following initial conditions, namely, the TLI maneuver, the initial state and the time of flight. A good initial estimate is essential for fast convergence to the final solution in the numerical generation of lunar free-return trajectories. Many analytical techniques are used as quick mission design and analytical tools for the initial estimate. These analytical techniques produce approximate closed-form solutions under certain assumptions. Typically, the actual Earth-moon system is approximated with point conic model or Lambert’s conic model [ 3 ], the matched conic model [ 4 ], the circular restricted three-body model (CRTBM) [ 5 ], and the pseudostate theory [ 6 ]. The initial estimate can be generated under the two-body Earth-spacecraft gravitational model by solving the Lambert’s problem [7]. Gibson determined that the matched conic

2

model can well approximate the results of precision trajectories and form a basis for an iteration scheme to converge upon the initial estimate for a more precise model in Ref. [8]. Li constructed a two-segment lunar free-return trajectory that retains the advantages of hybrid returns with the matched conic model in Ref. [9]. Based on the CRTBM, Schwaniger investigated the characteristics of the cislunar and circumlunar freereturn trajectories in Ref. [10 ]. Jesick and Ocampo developed a generation algorithm to construct a symmetrical single-flyby lunar free-return trajectory with matched conic and CRTBMs in Ref. [7]. Based on the pseudostate theory, Ramanan proposed an impact algorithm and a nonimpact algorithm for lunar transfer trajectories in Refs. [11,12]. Zhang developed an improved pseudostate technique to construct an accurate initial estimate for the lunar free-return trajectory in Ref. [13]. An iteration scheme is usually used to obtain a precision transfer trajectory by searching the final solution in a high-fidelity model based on the initial estimate. Green and Lewin proposed a gradient method for obtaining lunar free-return trajectories that consider constraints through a penalty function in Ref. [14]. The differential-correction method is one of the most common and efficient iteration schemes. Li [9] and Jesick [7] used this method and its variations in designing lunar free-return trajectories. Luo utilized the second-order state transition matrix technique to enlarge the convergence domain and to improve the convergence rate further in Ref. [15]. After converting this problem to a nonlinear programming, Peng used sequential quartic programming to design precision lunar free-return trajectories in Ref. [16]. Some global optimization algorithms have also been used for the search process. Ramanan proposed a numerical search technique through genetic algorithms with adaptive bounds to design lunar gravity assist transfers to geostationary orbits in Ref. [17]. In the case of a lunar free-return trajectory, the differential-correction method is generally used as the iteration scheme to search the final solution. Difficulties can arise in using this classic method. Given that differential equations are unstable or very stiff, they are highly sensitive to the initial conditions. Consequently, a good initial estimate is required to activate the iteration process. Otherwise, the forward integration can overflow, or the Jacobian matrix used in computing the corrections can be severely ill conditioned. However, a good initial estimate is often difficult to obtain through easy procedures without prior knowledge of the solution geometry. To solve the aforementioned problems, a numerical search method based on UKF parameter estimation (UKF-PE) is proposed in the present study. The algorithm effectively evaluates the Jacobian matrix precisely through its sigma points. Thus, without the need to perform any analytic differentiation, the calculation of the Jacobian matrix and other derivative terms that are often computationally difficult and that constitute the sources of numerical ill conditioning is omitted. Moreover, the difficulty of guessing the initial conditions decreases because the algorithm has large convergence ability. The algorithm allows the initial estimate to be generated under the two-body Earth-spacecraft model by solving the corresponding quasi-Lambert’s problem. Based on this initial estimate, the algorithm can easily converge to the final solution under a high-fidelity model in a few iterations by converting the original design problem to a new parameter estimation problem. Thus, the algorithm is employed to accurately and efficiently solve the problem of generating the lunar free-return trajectory.

STATEMENT OF THE PROBLEM Gravitational Model The initial estimate can be generated under the two-body Earth-spacecraft model, where the terminal attraction of the moon on the spacecraft is ignored, and the moon moves in a circle orbit around the Earth. Thus, the original generation of the lunar free-return trajectory is converted to a rendezvous or intersection problem with the spacecraft and the moon in the gravitational field of the Earth. Based on the initial estimate, a high-fidelity gravitational model is used to search the final solution. In addition to the central gravitational forces of the Earth and moon, the J 2 non-spherical perturbation of the Earth, the three-body perturbation of the Sun, and solar radiation pressure (SRP) are considered. The real positions of the moon and the Sun are obtained from the Jet Propulsion Laboratory DE405 ephemeris. The equations of motion in the J2000 ECI frame are described as

3

x = vx , y = v y , z = vz

(1)

2 3 § RE · § z 2 · º μ M ( x − x M ) μ S ( x − xS ) vx = − 3 «1 + J 2 ¨ − ¸ ¨1 − 5 2 ¸ » − rE « 2 © rE ¹ © rE ¹ » rM3 rS3 ¬ ¼ μ x L x − xS μ x − M3 M − S3 S + KCR Am s 4π cl rS3 rEM rES

(2)

μE x ª

3 § RE · § z 2 · º μ M ( y − yM ) μ S ( y − y S ) « 1 1 5 J + − − ¨ ¸» − ¨ ¸ 2 rE3 « 2 © rE ¹ © rE2 ¹ » rM3 rS3 ¬ ¼ μ y L y − yS μ y − M3 M − S3 S + KCR Am s 4π cl rS3 rEM rES

v y = −

μE y ª

μE z ª

2

2

3 § RE · § z2 « 1 1 5 J + − ¸ ¨ 2¨ rE3 « 2 © rE ¹ © rE2 ¬ μ z L μ z − M3 M − S3 S + KCR Am s 4π cl rEM rES

vz = −

· º μ M ( z − zM ) μ S ( z − zS ) − ¸» − rM3 rS3 ¹ »¼ z − zS rS3

(3)

(4)

where °­0, if rS ⋅ r < 0 and rS × r < rS RE K =® °¯1, else

(5)

In this study, the cylindrical shadow model is used for SRP. CR and Am are chosen to be 1.95 and 20 m 2 , respectively.

Constraint Analysis For the lunar free-return trajectory, if the injection time and transfer time are specified, the whole trajectory is determined by the initial state of the outbound trajectory. To meet the demands of human lunar missions, the constraints are imposed at three points of particular interest in the complete trajectory: the injection point, the perilune, and the reentry point. At the injection point, io and Ro are determined by the launch vehicle and the location of the launch site. Moreover, the injection point is chosen to be the perigee of the outbound trajectory that means γ o = 90 deg. At perilune, rp and tM are specified to guarantee favorable observation conditions and to complete the perilune brake. At the reentry point, γ re is always a negative value that is typically located at [ −7.5, −5.5] deg, similar to the Apollo reentry conditions. This angle should be selected such that the reentry into the Earth’s atmosphere is within the entry corridor. Moreover, hre is usually 120 km according to the altitude of the dense atmosphere.

PRELIMINARY DESIGN ALGORITHM The motion of a spacecraft in the Earth-moon system can be approximated with a two-body Earthspacecraft model. In this model, the trajectory of the spacecraft is a conic curve, which is referred to as a Keplerian orbit. The TPBVP of the Keplerian orbit has been studied intensively. In general, as shown in Figure 1, the Keplerian orbit boundary problem can be described with variables r1 , r2 , θ , γ 1 , γ 2 , and t12 . Let O be the center of the gravity field, and P1 , P2 be the initial and final points on the orbit, respectively. t12 is the flight time from P1 to P2 . r1 and r2 are the radial distances from P1 and P2 to O , respectively. θ is the transfer angle between the vectors r1 and r2 with respect to O . γ 1 and γ 2 are the flight-path angles at P1 and P2 , respectively.

4

Lambert’s problem is one of the most famous boundary problems for the Keplerian orbit, which has been successfully applied in the design of near-Earth and deep space orbits. However, in the design of the lunar transfer orbit, the problem of determining a Keplerian orbit based on r1 , r2 , γ 1 , and t12 is frequently encountered. This kind of problem is referred to as quasi-Lambert’s problem with fixed flight-path angle constraint. It determines a Keplerian orbit which connects two points that are assigned only with radius distances in a specified transfer time with a fixed flight-path angle at the initial point. An efficient method for solving this problem is proposed in Ref. [18]. The method is used to construct the initial estimate for free-return outbound trajectories in this study. Z

v2

γ2 P2

r2

O

θ

Y

t12

X

r1 P1

γ1

i

v1

Figure 1. Geometry of TPBVP in Keplerian orbit. Based on the aforementioned quasi-Lambert’s problem, a convenient preliminary design algorithm without a need to calculate gradients is proposed to construct an initial estimate for the lunar free-return outbound trajectory in the two-body Earth-spacecraft model. The algorithm is divided into three major steps: 1) Based on the preceding constraint analysis, r0 , γ 0 , TM , and tEM are specified. Moreover, the distance from the Earth to the moon rEM at TM is obtained from the DE405 ephemeris. Thus, the problem is converted to a quasi-Lambert’s problem with fixed flight-path angle constraint. 2) A free-return outbound trajectory that satisfies the imposed constraints is generated with the method in Ref. [18] for solving the quasi-Lambert’s problem proposed. Two solutions, which correspond to the twofold ambiguity in the orientation of the outbound orbital planes, can be generated. The initial state is transferred to the J2000 ECI frame by considering i0 . Thus, the three-dimension initial state x0 is generated. 3) Kepler’s equations are adopted to convert the spacecraft state x0 from the position and velocity to the orbital elements, and the initial estimate w0 = [vTLI, 0 , Ω TLI, 0 , ωTLI, 0 ]T of the target lunar free-return trajectory can then be constructed.

NUMERICAL TARGETING FOR THE FINAL SOLUTION Based on the initial estimate w0 , the problem of finding an appropriate numerical iteration scheme to T * * º¼ that corresponds to the high-fidelity gravitational model is search the final solution w * = ª¬vTLI Ω*TLI ωTLI studied in this section. For this problem, the free parameters are chosen to be

ªvTLI º w = ««Ω TLI »» «¬ωTLI »¼ 3×1 and the constraints to be satisfied are chosen to be

5

(6 )

ª( rp' − rp', t ) / rp', t º « » d = «( γ p' − γ p', t ) / γ p', t » « » «( γ re − γ ret ) / γ ret » ¬ ¼

(7 ) 3×1

Next, d is considered to be the output of a nonlinear transformation of w as d = F (w )

(8)

where F is the mapping function that in fact integrates numerically from the initial state to the final state. The problem is essentially solving a system of three nonlinear equations with three unknowns, F ( w ) = 0 .It can be stated as to find a vector w * such that F ( w * ) = 0 . This solution can be obtained as a fixed point of some function G by means of fixed point iteration

wk +1 = G ( wk ), k = 0,1,...

(9)

wk +1 = wk − J k-1d k , k = 0,1,...

(10)

A common and efficient method for solving the system of nonlinear equations, F ( w ) = 0 , is the differential-correction method. The generic iteration of the differential-correction method is where J k is the Jacobian matrix of the function F ( w ) evaluated at the kth iteration wk , i.e., T

T T T (11) J k = ª∇F1 ( wk ) ∇F2 ( wk ) " ∇Fn ( wk ) º ¬ ¼ This method is widely used for its fast convergence rate and accuracy. However, it may fail to converge if the initial estimate is not close enough to the final solution due to its property of having local convergence. Moreover, the common way to construct J k is to use the finite difference technique that may lead to singular or ill-conditioned matrix and hence make it more difficult to solve the original problem.

The UKF-PE method is proposed as an alternative way to solve this problem accurately and efficiently. Different from the differential-correction method in nature, let w be a 3-dimensional random variable with the mean wk and given covariance Pwk . A set of weighted points S = {Wi , Ȥ i } (sigma points) are properly chosen to make their mean wk and covariance Pwk by some sampling strategy. Each sigma point is propagated through the nonlinear transformation F to yield the set of projected sigma points L = {Wi , d i } . The estimated mean d k and covariance Pdk d k of the constraints can be approximated using a weighted sample mean and covariance of the projected sigma points, referring to extensive description of unscented transformation in Ref. [19]. Then a sequence {wk } can be constructed using the updating step of the UKF process wk +1 = wk − K k d k , k = 0,1,...

(12)

where K k is the Kalman gain matrix, , referring to extensive description of UKF-PE in Ref. [20]. Through updating wk step by step, the sequence {d k } approaches zero. When d k and its covariance Pd k d k is close enough to zero, the constraints imposed can be considered as satisfied, and hence wk is considered as the final solution w * . The whole process is illustrated in Figure 2. Though similar to the differential-correction method in form, the UKF-PE method is essentially based on the theory of probability. In addition to its deceptively simplicity, the UKF-PE method has two main advantages. First, the UKF-PE effectively evaluates the Jacobian precisely through its sigma point propagation, without the need to perform any analytic derivation or numerical differentiation. Therefore, the UKF-PE can be used to handle the cases with discontinuities or with non-analytical Jacobians. Second, the fact that the convergence ability of UKF-PE outperformed the differential-correction method makes the initial guessing a much simple task. To utilize the UKF-PE method, the original problem is converted to a new state-space representation, wk +1 = wk + rk 0 = F ( w k ) + ek

6

(13)

Pd1d1

Pw0

Pw2

d1

w0

w2 w*

Pd0 d0

d = F (w )

0

d2

d0

w1

Pd2 d2

Pw1

Figure 2. UKF parameter estimation. where w is the estimated parameter corresponding to a stationary process with the identity state transition matrix, d is the output that is expected to be zero, and F is the nonlinear measurement on wk . The process noise rk ~ N (0, R r ) drives the dynamic system, and the observation process is given by ek ~ N (0, R e ) . Note that rk and ek should be zero for the deterministic problem. However, for the sake of numerical stability, R r and R e are usually set to small values. Then, the UKF-PE method can be applied directly as an efficient method for learning the parameters, as illustrated in Figure 3. The scaling factor α determines the spread of the sigma points and is usually set to [10−4 , 1] . Moreover, β is used to incorporate prior knowledge of the distribution of w and usually set to zero. Two common options exist on how to choose Rkr , where η ∈ [ 0,1] , ρ RLS ∈ (0,1] . Note that stability can be improved through a proper scaling of the problem. In all the cases studied, the free parameters and constraints are scaled to achieve a value on the order of one to improve the stability of numerical searching. The computational complexity of UKF-PE is highly correlated with the total number of sigma points. This observation is uniquely determined by the sampling strategy while the dimension of the parameters is selected. In order to obtain minor computational complexity and robust convergence performances, the spherical simplex sampling [21] is chosen to be the sampling strategy in this study. Moreover, given that this problem is highly nonlinear, the scaled sampling [22] is combined with the spherical simplex sampling. Thus, the generation of sigma points is as follows: ­°W0 / α 2 + 1 − 1 / α 2 , i = 0 Wi ( m ) = ® 2 °¯(1 − W0 ) / α (n + 1), i = 1, 2,! n + 1 (m) 2 °­W0 + W0 + 1 + β − α , i = 0 Wi ( c ) = ® ( m ) °¯Wi , i = 1, 2,! n + 1

Ȥ 10 = [0], Ȥ11 = [−1 / 2W1( m ) ], Ȥ 12 = [1 / 2W1( m ) ]

­ j −1 °ª Ȥ0 º , i = 0 ° «¬0 »¼ ° j −1 °ª º ° Ȥi j » , i = 1,! , j Ȥi = ® « ( m) ° «¬ −1 / j ( j + 1)W1 »¼ ° º ° ª0 j −1 » ,i = j +1 °« (m) °¯ ¬« j / j ( j + 1)W1 ¼»

for j =2,! ,n

where W0 is the weight on the first point of spherical simplex sampling.

7

(14)

(15) (16)

(17)

w0 Pw0 R0r k =1

R0e

wk− = wk −1

Pw−k = Pwk −1 + Rkr −1 Wk |k −1 = wk− + α Pw−k Ȥ kn Dk |k −1 = H (Wk |k −1 )

d k = H ( wk− )

k = k +1 n +1

Pd k d k = ¦Wi ( c ) [ Di , k |k −1 − d k ][ Di , k |k −1 − d k ]T + Rke−1 i =0 n +1

Pwk d k = ¦Wi ( c ) [Wi , k |k −1 − wk ][ Di , k | k −1 − d k ]T i=0

K k = Pwk dk Pd−k1dk wk = wk− − K k d k Pwk = Pw−k − K k Pdk d k K kT

option 1: Rkr = (1 − η ) Rkr −1 + η K k H ( wk ) H ( wk ) K kT T

−1 option 2 : Rkr = ( ρ RLS − 1) Pwk

Rke = Rke −1

Figure 3. UKF-PE flowchart.

NUMERICAL EXAMPLES Some simulations are carried out in the high-fidelity model to illustrate the accuracy and convergence performance of the free-return generation procedure. A comparison is made with trajectories that satisfy the same constraints between two iteration schemes: the differential-correction method and the UKF-PE method. Note that the differential-correction method in this study uses the Broyden technique that is free of calculating the Jacobian matrix in Ref. [23]. Lunar free-return trajectories include circumlunar and cislunar trajectories and departure form posigrade and retrograde Earth orbits. The analysis presented in this paper is limited to the circumlunar free-return trajectory, for which perilune occurs at the far side and passed all the way around the moon. Since the generation procedure is greatly influenced by the Earth-moon geometry, four special cases in a month are chosen, as listed in Table 1. To begin with, the initial estimates of these four cases under the two-body Earth-spacecraft model can be generated by solving the corresponding quasi-Lambert’s problem, as shown in Table 2. We can observe that the initial estimate generated through the simple and mature preliminary design algorithm is relatively poor compared with the final solution. The maximum error on the injection velocity is approximately 3 m/s, and the maximum error on the angles error is approximately 2 deg. The initial estimate increases the difficulty to converge for the search algorithm because of high sensitivity of this problem. In fact, for each

8

case, two solutions are obtained corresponding to departure from posigrade and retrograde Earth orbits. Without loss of generality, the retrograde trajectories are only discussed in this study.

Table 1. Initial conditions Case TTLI , Jd

1

2

3

4

2455353

2455360

2455366

2455373

rEM , km

~402000

~392000

~372000

~370000

t EM , h

71.0

60.0

60.0

65.0

ho , km

263.5

263.5

263.5

263.5

io , deg

24.0

30.0

24.0

35.0

hp , km

120.0

120.0

120.0

120.0

hre , km

120.0

120.0

120.0

120.0

γ re , deg Remarks

-6.0 Apogee of the moon’s orbit

-6.0 Smallest lunar declination

-6.0 Largest lunar declination

-6.0 Perigee of the moon’s orbit

Based on the initial estimate, the differential-correction method and the UKF-PE method are used to search the converged final solution. Numerical integration is performed in the high-fidelity model with a Runge-Kutta-Fehlberg 12 (13) variable-step integrator. The stopping criterion used in single performance of the program is chosen to be 1/ 2

§ 3 · f ( d ) = ¨ ¦ di2 ¸ © i =1 ¹

<ε

(18)

where the error tolerance ε = 1e − 5 . Therefore, we check that iterates converge to a limit and that this limit is a solution of the problem. Moreover, the maximum iteration of each method is set to 300. As for the detailed settings of UKF-PE, the covariance Pw0 is set to I , R0r is set to 1e − 4 I , and R0e is set to 1e − 20 I . The weight on the first sigma point W0 is set to 0.5. The scaling factor α is set to 5e − 4 , and β is set to 0. The updating method of R0r is chosen to be option 1 in Figure 3, and the weighting factor η is set to 0.5. Performances of the differential-correction method and the UKF-PE method are explored numerically. The differential-correction method fails to obtain a converged solution for the four cases above. The UKFPE method is able to find the converged final solution in a few iterations and obtains robust convergence performances. The initial estimate (IE), final solution (FS), value of constraint function ( f (d ) ), and number of iterations needed to converge to the solution (Iter) are analyzed, as shown in Table 2, and the errors of targeting parameters are shown in Table 3.

Table 2. Targeting performance of UKF-PE Case 1, IE Case 1, FS Case 2, IE Case 2, FS Case 3, IE Case 3, FS Case 4, IE Case 4, FS

vTLI , m/s

ȍ TLI , deg

ωTLI , deg

10896.91 10894.42 10902.74 10899.75 10901.12 10899.09 10905.10 10902.04

347.746 347.488 27.012 28.375 166.280 168.922 303.206 304.474

190.835 192.225 248.028 247.735 191.062 189.140 142.659 142.534

9

f (d ) 52.839 4.32e-6 41.946 5.30e-7 51.003 4.90e-7 44.935 1.47e-6

Iter 18 11 21 12

Table 3. Errors of the targeting parameters of UKF-PE

Case 1, FS Case 2, FS Case 3, FS Case 4, FS

rp' , m

γ 'p , deg

rre , m

7.41 -0.44 -0.38 2.62

4e-5 1e-5 1e-5 -1e-5

10.95 3.12 2.87 -0.41

We can conclude from Table 2 that the UKF-PE method converges to the final solution in less than 30 iterations and shows good convergence behavior. Table 3 shows that the error on the achieved altitude is less than 15 m, and the error on the achieved flight-path angle is less than 1e-4 deg, which means the constraints can be satisfied accurately to guarantee effectiveness of the simulation results. Therefore, the UKF-PE method decreases the difficultly of guessing a good approximation encountered in the differentialcorrection method. Moreover, the method can be utilized as an alternative tool to generate the lunar freereturn trajectory robustly and efficiently. The feasible and converged lunar free-return trajectories via the algorithm proposed previously are shown in Figure 4.

Figure 4. Lunar free-return trajectory in the J2000 ECI frame. Table 4. Convergence domains in the one-dimensional search space Method Disturbance

Differential-correction ΔvTLI ,m/s Δȍ TLI , deg

Case 1 Case 2 Case 3 Case 4 Maximum Minimum Average

13.0 9.5 11.1 19.3 19.3 9.5 13.23

6.1 1.0 3.3 3.4 6.1 1.0 3.45

ΔωTLI , deg

UKF-PE ΔvTLI , m/s

Δȍ TLI , deg

ΔωTLI , deg

2.1 3.5 3.5 1.3 3.5 1.3 2.60

35.7 39.6 35.5 38.9 39.6 35.5 37.43

14.1 4.4 5.1 3.4 14.1 3.4 6.75

5.4 7.2 7.3 8.3 8.3 5.4 7.05

To examine the convergence domain of the UKF-PE method, a disturbance is added on one coordinate component of the final solution, whereas the other two components remain invariant. From the disturbed point, these two methods are used to search the original final solution. Let the magnitude of the disturbance be continuously increased until the searching process fails to obtain a converged solution. Then, the

10

convergence domains of these two methods in the one-dimensional search space for the disturbed component are obtained. The statistics of the convergence domains are shown in Table 4. According to the preceding statistical data, the geometric structure of the convergence domain has no specific regular pattern because of the high nonlinearity of the lunar free-return trajectory, and the convergence domains of the components are greatly influenced by the Earth-moon geometry. Moreover, the UKF-PE method obtains a large convergence domain for these examples. In the worst situation, the convergence domain of the differential-correction method for ȍ TLI and ωTLI is approximately 1 deg, which means a very good initial guess is necessary. Note that the UKF-PE method does not require more iterations than those of differential-correction method. To further validate the convergence ability of the UKF-PE method, a bounded disturbance is determinately added on each component of the final solution, where ΔvTLI ∈ [ 0,35] , Δȍ TLI ∈ [ 0,14] , ΔωTLI ∈ [ 0,5] . Let I v , I Ω , and Iω be the number of points in which the domains of ΔvTLI , Δȍ TLI , and ΔωTLI defined earlier are discretized, respectively. The three-dimensional search space has been sampled with I v = 18 , I Ω = 14 , and Iω = 10 . Thus, in total 2520 disturbed first guess orbits are generated. Based on these first guesses, these two methods are used to search the final solution. Without loss of generality, the analysis in this study is limited to case 1. The convergence performances are shown in Figure 5, in which ' ' the converged points are plotted, and ΔvTLI ∈ [ 0,1] are normalized units. In ∈ [ 0,1] , Δȍ 'TLI ∈ [ 0,1] , Δω TLI order to represent the distribution of points, the percentages of the converged points at different radical ' ' ' distances are calculated. The radical distance of each disturbed point x ( ΔvTLI , ΔΩTLI , ΔωTLI ) is defined as

( Δv ) + ( ΔΩ ) + ( Δω ) ≤ r , converged} , N = { y | ρ ≤ r } , ρx =

2 ' TLI

' 2

' 2

(19)

where ri +1 = ri + 0.1 , r1 = 0.2 , Consider collections M i = { x | ρ x i i y i i = 1, 2, !,16 . Let mi , ni be the numbers of elements in the collections M i , N i , respectively. Therefore, the percentage of the converged points is defined as Pi = mi / ni , shown in Figure 6. From Figure 5, we can conclude that the differential-correction method converges in 159 (6.3%) cases, and the UKF-PE method converges in 1968 (78.1%) cases. Therefore the convergence domain of the UKFPE method is, on the average level, approximately 12.4 times as large as that of the differential-correction method for these examples. From Figure 6, as the radical distance increases, the percentage of the converged points for the differential-correction method significantly decreases, while the percentage for UKF-PE remains relatively high. Note that when the disturbed points lie within the unit sphere ρ x ≤ 1 , the percentage of convergence for UKF-PE reaches 93.1%. The simulation results also demonstrate that the disturbed points that fail to converge with UKF-PE are mostly located at a relatively far distance ρ x > 1 exceeding the convergence ability of this method with the selected values of its adjustable parameters.

Differential-correction

UPE

ΔωTLI,deg

1

0.5

'

0.5

'

ΔωTLI,deg

1

0 1

0 1

1 0.5

ΔΩ 'TLI,deg

0.5 0 0

1 0.5

Δv'TLI, m/s

ΔΩ 'TLI,deg

0.5 0 0

Δv'TLI, m/s

Figure 5. Convergence domains in the three-dimensional search space.

11

100%

Differential-corretiion UPE

90% 80% 70%

P

60% 50% 40% 30% 20% 10% 0 0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

ρ

Figure 6. Percentages of the converged points. To examine the convergence rate of the UKF-PE method, several fixed disturbances ΔvTLI are added on the first coordinate component of the final solution for case 1. From these disturbed points, these two methods are used to search the original final solution. The iterative processes of these two methods are shown in Figure 7. We can conclude that when the disturbance is small, the UKF-PE method converges more slowly than the differential-correction method. As the disturbance increases, the differentialcorrection method may fail to converge. The simulation results partly demonstrate the advantage of the differential-correction method with fast convergence under the condition that the initial estimate is close to the final solution. Therefore, we can conclude that these two types of methods can be combined for solving the problem resulting in a good balance between convergence domain and rate. Differential-correction

2

UPE

2

10

10

ΔvTLI=1

ΔvTLI=1

ΔvTLI=5

0

10

ΔvTLI=5

1

10

ΔvTLI=10

ΔvTLI=10

ΔvTLI=15

ΔvTLI=15

0

10

-2

10

-1

-4

f(p)

f(p)

10

10

-2

10

-3

10

-6

10

-4

10 -8

10

-5

10

-10

10

-6

0

5

10 Iterations

15

10

20

0

5

10 Iterations

15

20

Figure 7. Iterative processes of methods.

CONCLUSION This paper demonstrates that the UKF-PE technique represents a convenient, efficient and accurate method for designing lunar free-return trajectories. The method is intuitive and straightforward; it is without the need of the derivation of the gradient matrixes; it only requires the initial guess generated under a two-body Earth-spacecraft model. Thus, the difficulty of the design of lunar free-return trajectories is greatly reduced. The lunar free-return trajectory generated can well satisfy the constraints at injection, perilune, and Earth-entry conditions encountered in lunar exploration missions. The outstanding

12

performances of this method have been verified through extensive numerical simulations. Therefore, the method proposed provides an alternative way to design the lunar free-return trajectory conveniently and efficiently.

1 Orloff, R. W., “Apollo by the Numbers: A Statistical Reference,” U.S. Government Printing Office, Washington, D.C., 2000, pp. 34, 36-37, 74-75,92-93,115-116, 163,187,215,243. 2

Low, G. M., “Apollo 11 Mission Report,” NASA MSC-00171, Nov. 1969.

3

Battin, R. H., and Vaughan, R. M., “An Elegant Lambert algorithm,” Journal of Guidance, Control, and Dynamics, Vol. 7, No. 6, 1984, pp. 662-670. 4

Egorov, V. A., “Certain Problems of Moon Flight Dynamics,” Russian Literature of Satellites, Pt. 1, Vol. IED, International Physical Index, New York, 1958, pp. 107-175. 5

Bruno, A. D., The Restricted 3-body Problem: plane periodic orbits, Walter de Gruyter, New York, 1994, pp. 95-107.

6

Wilson, S. W., “A Pseudostate Theory for the Approximation of Three-Body Trajectories,” TRW Note No. 69-FMT765 (11176-H304-R0-00), 1969. 7

Jesick, M., and Ocampo, C., “Automated Generation of Symmetric Lunar Free-Return Trajectories,” Journal of Guidance, Control, and Dynamics, Vol. 34, No. 1, 2011, pp. 98-106. doi: 10.2514/1.50550

8

Gibson, T., “Application of the Matched Conic Model in the Study of Circumlunar Trajectories”, NASA Project Apollo Working Paper 1066, Feb. 1963.

9

Li, J., and Baoyin, H., “Analysis of Two-Segment Lunar Free-Return Trajectories,” Journal of Spacecraft and Rockets, Vol. 52, Special Section on Numerical Simulation of Hypersonic Flows, 2015, pp. 183-195. doi: 10.2514/1.A32874

10

Schwaniger, A. J., “Trajectories in the Earth-Moon Space with Symmetrical Free Return Properties,” NASA TN D1833, June 1963. 11

Ramanan, R. V., “Intergrated Algorithm for Lunar Transfer Trajectories Using a Pseudostate Technique,” Journal of Guidance, Control, and Dynamics, Vol. 25, No. 5, 2002, pp. 946-952. doi: 10.2514/2.4968 12

Ramanan, R. V., and Adimurthy, V., “Nonimpact Lunar Transfer Trajectories Using the Pseudostate Technique,” Journal of Guidance, Control, and Dynamics, Vol. 28, No. 2, 2005, pp. 217-225. doi: 10.2514/1.7607 13

Zhang, H., Luo, Q., and Han, C., “Accurate and Fast Design Algorithm for Free-Return Lunar Flyby Trajectories,” Acta Astronautica, Vol. 102, 2014, pp. 14-26. doi:10.1016/j.actaastro.2014.05.015 14 Green, B. S., and Lewin, N., “A Gradient Method for Obtaining Circumlunar Trajectories”, AIAA Paper 63-401, Aug. 1963. 15

Luo, Q., Yin, J., and Han, C., “Design of Earth-Moon Free-Return Trajectories,” Journal of Guidance, Control, and Dynamics, Vol. 36, No. 1, 2013, pp.263-271. doi: 10.2514/1.55910 16

Peng, Q., Shen, H., Li, H., “Free Return Orbits Design and Characteristics Analysis for Manned Lunar Mission,” Science China Technological Sciences, Vol. 54, No. 12, 2011, pp. 3243-3250. doi:10.1007/s11431-011-4622-7 17 Ramanan, R. V., and Adimurthy, V., “Precise Lunar Gravity Assist Transfers to Geostationary Orbits,” Journal of Guidance, Control, and Dynamics, Vol. 29, No. 2, 2006, pp. 500-502.

doi: 10.2514/117469

13

18

Luo, Q., Meng, Z., and Han, C., “Solution Algorithm to a Quasi-Lambert’s Problem with Fixed Flight-Direction Angle Contraint”, Celestial Mechanics and Dynamical Astronomy, Vol. 109, No.4, 2011, pp. 409-427. doi: 10.1007/s10569-011-9335-5 19

Julier, S. J., and Uhlmann, J. K., “Unscented filtering and nonlinear estimation”, Proceedings of the IEEE, IEEE, Vol.92, Mar 2004, pp. 401-422. doi: 10.1109/JPROC.2003.823141 20 Wan, E. A., and R. van der Merwe, “The Unscented Kalman Filter,” Kalman Filtering and Neural Networks, Wiley, New York, 2001, pp. 221-280 21

Julier, S. J., “The Spherical Simplex Unscented Transformation,” American Control Conference, 2003. Proceedings of the 2003, IEEE, Vol.3, June 2003, pp. 2430-2434. doi: 10.1109/ACC.2003.1243439

22

Julier, S. J., “The Scaled Unscented Transformation,” American Control Conference, 2002. Proceedings of the 2002. IEEE, Vol. 6, May 2002, pp. 4555-4559. doi: 10.1109/ACC.2002.1025369 23 Broyden, C. G., "A Class of Methods for Solving Nonlinear Simultaneous Equations," Mathematics of computation, Vol. 19, 1965, pp. 577-593

14