# Low-Thrust Trajectory Design Using Finite Fourier Series Approximation of Pseudoequinoctial Elements

Authors: Wanmeng Zhou; Haiyang Li; Hua Wang; Ran Ding
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

A new low-thrust trajectory design method is proposed that is based on the finite Fourier series method with pseudoequinoctial elements rather than the more common cylindrical coordinate components. The bijection relation between the elements and control variables is ensured by introducing an additional equality constraint derived from the angular momentum conservation. The guidance law and on-line control variables are obtained by applying inverse dynamics and the framework of inverse simulation technology, respectively. The pseudoequinoctial finite Fourier series method has the advantages of both the Fourier series and the perturbation analysis methods. For two-body problems, three cases were studied: the Earth to Mars, 1989ML, and Tempel-1 missions. Regarding the design of a rendezvous trajectory with a large inclination angle and a high eccentricity rate, this method yields a broader range of feasible results than the traditional Fourier series method. The circular restricted three-body problem was solved for the first time using the pseudoequinoctial finite Fourier series method combined with the patched conics method. The low-thrust Earth-Moon transfer was analyzed, and the results show that this method improves window analysis efficiency and guarantees precision of the initial geocentric trajectory for the low-thrust transfer.


---

## Body

## 1. Introduction

In recent years, low-thrust spacecraft have increasingly attracted interest in the field of engineering because of their long-lasting and fuel-saving features. These include the National Aeronautics and Space Administration’s (NASA) Deep Space 1 [1] and Dawn [2], the European Space Agency’s (ESA) SMART-1 [3], and BepiColombo [4], a joint mission between ESA and Japan Aerospace Exploration Agency (JAXA). The trajectory design of these low-thrust vehicles has also attracted much attention. In a global trajectory optimization competition (GTOC 2), the trajectory design problem for low-thrust vehicles was proposed to explore asteroids and issues associated with them [5]. The challenges in the search for the optimal trajectory of a low-thrust, deep-space exploration vehicle are greater than those of an impulsive-thrust design due to features such as the long time of flight (TOF) and multiple solutions of the feasible flight trajectory. Given the initial and final position, velocity, and TOF, the thrust profile can be determined by optimal fuel or energy under maximum thrust amplitude, which is also called a two-point boundary value (TPBV) problem [6, 7].There are direct and indirect methods for solving TPBV problems. In indirect methods [8], according to Pontryagin’s maximum principle, the derivative of the Hamiltonian function with respect to the control input equals zero. The control input can then be represented by both state and costate variables [9]. Substituting the control back into the differential equations of the costate and state variables, the simultaneous equations are then solved based on the boundary and transversal conditions. The disadvantage of these methods is that it is difficult to assign the costates’ appropriate initial values when solving the problem, for their unapparent physical significances. In direct methods, such as point collocation [10] and the pseudospectral method [11], either the control variables or both the state and control variables are discretized to form a set of nonlinear equations with unknown parameters that are solved by nonlinear programming algorithms. Although these methods do not require the introduction of costate variables, the discrete points may become increasingly dense during the iteration process, which makes the calculation time consuming and potentially unstable. Other variants of direct methods are the shape-based methods, in which the flight states are suitably represented by parameterized nonlinear functions with certain assumptions and use inverse dynamics to retrieve the control corresponding to the assumed functional form. These methods require fewer collocation points and design parameters to obtain higher computational efficiency. Based on previous discussions of the optimal direction of a fixed thrust [12–14], early shape-based methods on two-dimensional problems were generally based on the tangential thrust assumption and the analytic results were typically expressed in polar coordinates [15–17] such as in the logarithmic spiral method, Forbes spiral method, exponential sinusoid method, and inverse polynomial method (IP). As the number of shape parameters increases, the number of constraint conditions that can be satisfied also increases. The design parameters can be derived directly from the constraint conditions without numerical iterations. To solve three-dimensional (3D) problems, there are several methods such as IP in cylindrical coordinates [18], the pseudoequinoctial method (PE) [19] and the spherical shaping method [20, 21] on the basis of the perturbation functions, and the hodograph method which considers the velocity variable as the parameterized shape [22]. The IP method cannot ensure the proper amplitude of the thrust, while PE and spherical shaping methods require additional design parameters to satisfy the TOF or thrust constraints [20]; neither can these methods provide the motion states directly related to the transfer time. For the hodograph method, an appropriate combination of basic functions should be selected first to obtain the optimal velocity.The finite Fourier series method (FFS) is a new method to overcome these deficiencies. It applies time as an independent variable and represents the position variables with an orthogonal polynomial with a finite number of trigonometric terms, i.e., Fourier series [23–25]. The orthogonal polynomials have been proven to solve optimal control problems successfully [26, 27]. The Fourier series in particular can represent any periodical function and could even theoretically represent the switching-type function with sufficient trigonometric function terms [28], which introduces more shape categories into the trajectory design. On the other hand, most research on the shape-based methods is applied on two-body problems, but few are applied for circular restricted three-body problems (CRTBPs). FFS combined with the simplified analysis method is the only shape-based method proposed so far to solve CRTBPs [29].Although FFS produces excellent results in transfer trajectory design between coplanar circular orbits [30], the increasing eccentricity and inclination angle of the transfer orbit will cause the total velocity increment to gradually diverge from the near-optimum solution and there will be more infeasible points in the launch window graph [31]. The first major reason is because there are a fixed number of terms selected which leads to a relatively fixed fitting shape. Another reason is because the expression of 3D-FFS in a cylindrical coordinate system has limited approximation capability. During the long-time transfer, the motion in the z-direction fluctuates periodically due to changing eccentricity, inclination angle, and the semimajor axis, which is difficult to match. In this study, the PE elements were introduced into the modified scheme of the FFS, denoted as PE-FFS. The PEs are similar to modified equinoctial elements (MEEs), which can provide smooth variations of instantaneous motion states compared to the z-axis components of a cylindrical coordinate system [32]. For two-body problems, different numbers of terms in each PE are first chosen for the different relative positions of the initial and final radius vectors. A feasible shape library expressed by combinations of these terms is then used to produce a wider range of feasible solutions. For CRTBPs, the new method proposed is the PE-FFS combined with the patched conics method. The capture situation and feasibility of the transfer trajectory are comprehensively analyzed through the state change at the capture point of the lunar sphere of influence (LSOI).This paper is organized as follows: First, the dynamic equations in the cylindrical coordinate system and the Gaussian perturbation equations of the modified equinoctial orbit elements are established separately under the two-body assumption. Then, the inverse dynamic model and inverse simulation technique are proposed to give the relationship between the state and control variables, and an additional condition is proposed to form the PE-FFS method. The third section presents the results of studies of the two-body rendezvous cases, which include transfers to Mars, an asteroid, and a comet. The results from the PE-FFS method are compared with those from the other methods. Finally, the CRTBP for the Earth-Moon transfer is investigated. The quick trajectory design and the launch window analysis for the Earth-Moon mission are realized by combining the PE-FFS and patched conics method.


## 2. Pseudoequinoctial Finite Fourier Series Method

### 2.1. Dynamic Modeling of a Two-Body System

A typical two-body problem [19] can be described as
(1)r¨+μr3r=f,where μ refers to the gravitational constant of the central body and f refers to the acceleration produced by the external force. In traditional shape-based methods, the state variable is represented by the polar radius r, the phase angle θ, and the height z along the z-axis. The components of equation (1) [32, 33] can be expressed as
(2)r¨−rθ̇+μrs3=fr,rθ¨+2ṙθ̇=fθ,z¨+μzs3=fz,where s=r2+z2 and fr,fθ,fz refer to the acceleration components in the local cylindrical coordinates. The definition of directions in the cylindrical coordinate system is shown in Figure 1.Figure 1
Diagram of the inertial and local cylindrical coordinate systems.The component in thez-axis periodically changes with the phase angle, and its peak and trough values are jointly determined by the semimajor axis, eccentricity, and inclination angle of the orbit. The schematic diagram of the component is shown in Figure 2.Figure 2
Schematic diagram of changes in thez-component during a multirevolution transfer.In the cylindrical coordinate system, if the inclination angle and eccentricity are small and there are few transfer cycles, an ordinary polynomial is suitable for the initial estimation of the component along thez-axis. As the inclination angle, eccentricity, and number of transfer cycles increase, convergence along the z-axis will become more difficult and the deviation in total velocity increment will increase gradually [31]. However, the orbit elements are slow variables to solve the problem. To avoid the singularity caused by zero eccentricity and zero inclination angle, the modified equinoctial elements are introduced and defined as follows [19]:
(3)p=a1−e2,f=ecosω+Ω,g=esinω+Ω,h=tani2cosΩ,k=tani2sinΩ,L=ω+Ω+θ,where a, e, i, ω, Ω, and θ denote the classical orbit elements. The two-body dynamic equation is expressed by the modified equinoctial elements in the Gaussian perturbation equation [19]:
(4)ṗ=2pp/μwf,ḟ=pμsinLfr+p/μw1+wcosL+fft−p/μwhsinL−kcosLgfn,ġ=−pμcosLfr+p/μw1+wsinL+gft+p/μwhsinL−kcosLffn,ḣ=s2pμcosL2wfn,k̇=s2pμsinL2wfn,L̇=p/μwhsinL−kcosLfn+μpw/p2,where s2=1+h2+k2, w=1+fcosL+gsinL, and fr,ft,fn denote the three acceleration components in the radial-transverse-normal coordinate system. The radial axis is the direction moving from the central gravitational body to the spacecraft, the normal axis is the direction of the trajectory angular momentum, and the transverse axis is determined by the right-hand rule. The components of the osculating orbit position rrx,ry,rz can be expressed by the modified equinoctial elements [19] as
(5)rx=rs2cosL+h2−k2cosL+2hksinL,ry=rs2sinL−h2−k2sinL+2hkcosL,rz=2rs2hsinL−kcosL,where r is the magnitude of the position vector, and r and s can be expressed as
(6)r=p1+fcosL+gsinL,s2=1+h2+k2.


### 2.2. Calculation of Guidance Law and On-Line Control Inputs

#### 2.2.1. Guidance Law on Inverse Dynamics

One critical step of the shape-based design method is to represent the control variables with states. This process can be expressed in the form of an inverse model [34] when a typical system is considered:
(7)ẋ=Fx,t+Gx,tut,y=Hx,t.The derivative of the second term in equation (7) can be given by
(8)ẏ=∇Hx,tẋ=LfHx+LgHxut,where LfHx and LgHx, also called Lie derivatives [35], can be expressed as
(9)LfHx=∂Hx∂x⋅Fx,LgHx=∂Hx∂x⋅Gx.WhenLgHx≠0, it leads to the control variable by getting a single derivative:
(10)ut=LgHx−1ẏ−LfHx.Otherwise, it requiresk derivatives until LgLfk−1Hx≠0, so that
(11)ut=LgLfk−1Hx−1yk−LfkHx,where LfkHx and LgLfk−1Hx are defined, respectively, as
(12)Lfk−1Hx=LfLfk−2Hx=∂Lfk−2Hx∂xFx,LgLfk−1Hx=LgLfk−1Hx=∂Lfk−1Hx∂xGx.For the components of the cylindrical coordinatesr,θ,z, the first-order and second-order derivatives can be calculated; for the orbit elements, ṙ and r¨ are derived from equation (5). Derivatives of both dynamic models are then substituted into equation (11) to obtain the control variable ut.


### 2.3. On-Line Control Based on Inverse Simulation

In a practical situation, the system is affected by external disturbances. A closed-loop control is required for the practical system. All the expected state variables are directly expressed with the independent time variable. Thus, the inverse simulation technique is quite suitable for the calculation of the practical control variables.The discrete system is(13)xtk+1=xtk+∫tktk+1ẋtk,uτdτ,ytk+1=Hxtk+1,where Hxtk+1 refers to the observation function. According to the current states obtained in the real-time system and the expected states, the expected control variable udtk can be obtained through iteration of the following [35]:
(14)un+1tk=untk−Jn−1Hxntk+1−ydtk+1,where J denotes the Jacobian matrix, which can be expressed as
(15)Jn=∂yerrorn∂u=yerrornun+δu−yerrornun−δu2δu,yerrornu=Hxntk+1−ydtk+1.The whole simulation framework is shown in Figure3.Figure 3
On-line control loop under the inverse simulation framework.The expected outputydtk+1 and the current state xtk are obtained based on the shaped trajectory and the measured value ytk, respectively. The expected control variable udtk is then calculated through the inverse simulation to form a closed loop.


### 2.4. Pseudoequinoctial Elements Based on the Finite Fourier Series

The state variable is represented in the form of a Fourier series:(16)λt=kλ2+∑n=1nλkaλncosnπTt+kbλnsinnπTt,where λ denotes the PE elements p~,f~,g~,h~,k~,L~ represented by the FFS, ka and kb represent the coefficients of each trigonometric term, nλ represents the number of the terms, and T denotes the TOF.


#### 2.4.1. Calculation of the Initial Parameters

Based on the boundary conditions, the first two coefficients of the sine and cosine terms in each element expression can be calculated as [31]
(17)ka1=λ0−λf2−∑n=3nλkann≥3,n∈odd,ka2=λ0+λf−kλ2−∑n=3nλkann≥4,n∈even,(18)kb1=Tλ̇0−λ̇f2π−∑n=4nλnkbnn≥3,n∈odd,kb2=Tλ̇0+λ̇f4π−12∑n=4nλnkbnn≥4,n∈even,where λ0,λf denote the initial and final PE elements, respectively, and λ̇0,λ̇f represent the derivatives of the initial and final elements, respectively. For smooth transfer processes, the elements and their change rate at the beginning and end of the transfer orbit should be consistent with those of the initial and final orbits.After imposing the boundary conditions, the number of remaining unknown variables in equation (16) is 2nλ−3. The initial value of the unknown coefficients can be obtained from a series of collocation points that are in accordance with the initial estimation curves. The typical initial estimate can be obtained by a third-order polynomial that satisfies the boundary conditions:
(19)λinitt=λ̇0T+λ̇fT+2λ0−λftT3+3λf−λ0−2λ̇0T−λ̇fTtT2+λ̇0t+λ0.In this third-order polynomial curve, discrete points such as the Legendre-Gauss points and Legendre-Gauss-Lobatto points are selected as the collocation points.ndp refers to the number of points. By combining equations (16) and (19), the nonlinear equation set can be expressed as
(20)λinitndp×1=Trigndp×2nλ−3Kλ2nλ−3×1+Boundndp×1,where Kλ denotes the unknown Fourier series coefficients, and Trig and Bound are expressed as
(21)Trig=12cos3πTt⋯cosnλπTtsin3πTt⋯sinnλπTtnppr×1,Bound=ka1cosπTt+ka2cos2πTt+kb1sinπTt+kb2sin2πTtnppr×1.Eventually, the Moore-Penrose pseudoinverse [35] is applied to equation (20), through which the unknown coefficients Kλ can be expressed as
(22)Kλ=TrigTTrig−1TrigTλinit−Bound.


#### 2.4.2. Practical Velocity and Acceleration Determination

After determining the initial values of all the coefficients, the next step is to calculate the thrust acceleration. According to the inverse simulation theory,ṙ and r¨ deduced from equation (5) can be transformed into the control acceleration u. The spacecraft velocity v=ṙ is expressed as
(23)v=∂r∂x∂x∂t+∂r∂L∂L∂t,where x refers to p,f,g,h,k. The osculating velocity vosc can be expressed as
(24)vosc=∂r∂LL̇osc=∂r∂Lμpr2.The actual velocity can also be expressed as(25)v=vosc+vgau,where vgau is the gauge velocity.By substituting equations (5) and (6) into equation (24), the osculating velocity component can be expressed as
(26)vx=−1p/μs2sinL+h2−k2sinL−2hkcosL+g−2fhk+h2−k2g,vy=−1p/μs2−cosL+h2−k2cosL+2hksinL−f+2ghk+h2−k2f,vz=2p/μs2hcosL+ksinL+fh+gk.When the dynamic equation satisfies the Gaussian perturbation equation, the osculating velocity and actual velocity coincide at each time, i.e., the osculating condition is satisfied andvgau≡0. The position and velocity can be obtained from equations (5) to (26), and vice versa. Thus, the modified equinoctial elements are equivalent to the osculating position and velocity. The PE elements, in the form of the Fourier series, are applied in order to ensure the bijection relationship between the motion states and pseudoelements. The position that satisfies equation (5) and the velocity at that time is
(27)v=∂r∂x~∂x~∂t+∂r∂L~∂L~∂t,where x~ refers to p~,f~,g~,h~,k~. In order to satisfy angular momentum conservation, an additional constraint is introduced and given by
(28)∂L~∂t=μp~r2.Thus, based on equations (5) and (27) and the hypothetical condition equation (28), the bijection relationship between the pseudoelements and contemporary motion states is obtained.In the calculation ofr¨, it is necessary to obtain a further time derivative of equation (24):
(29)r¨=∂vosc∂x~∂x~∂t+∂vosc∂L~∂L~∂t+∂vgau∂x~∂x~∂t+∂vgau∂L~∂L~∂t.For the modified equinoctial elements,vgau≡0, equation (29) will not provide any information on the external acceleration, i.e., r¨=−μ/r3r. However, it is impossible for pseudoelements to ensure that the gauge velocity is exactly zero. Equation (29) can be used to calculate the actual acceleration. The change in L~ with time is faster than in other elements. It can be deduced from equation (3) that when the external acceleration is lower than the local gravitational acceleration, the approximation condition is
(30)∂x~∂t≪∂L~∂t.As∂r/∂x~ and ∂r/∂L~ are of the same order, the pseudoelements should then satisfy
(31)∂r∂x~∂x~∂t≪∂r∂L~∂L~∂t.Then, the total velocity could be expressed asṙ≈∂r/∂L~∂L~/∂t and the osculating velocity and the position could be expressed as
(32)vosc=∂r∂LL̇osc≈∂r∂L~∂L~∂t,(33)r¨=∂v∂x~∂x~∂t+∂vgau∂L~L̇osc+∂vosc∂L~L̇osc.Thus, the osculating condition is basically satisfied. It can be inferred from equations (5), (27), and (33), according to the inverse simulation structure, that the thrust acceleration is expressed in the following form:
(34)F=Fλt,λ̇t,λ¨t.


#### 2.4.3. Pseudoequinoctial Fourier Series Method

Eventually, within the transfer timet0−tf, the constraint and velocity increment can be separately expressed as
(35)c1F2−Fmax≤0,c2∂L∂t−L̇osc=0,(36)Δv=∫t0tfF2dt.The rectangular quadrature rule is utilized in this paper. In addition, to ensure a prograde trajectory, an additional constraint for the orbit direction can be added as(37)−c3L̇osc≤0,where c1−c3 denote the relaxation factors. The relationship between the fuel cost and final velocity increment is expressed as [18]
(38)mf=m01−exp−ΔvIspg0,where Isp denotes the specific impulse and g0 denotes the standard gravitational acceleration at sea level. According to the constraints of equations (35) and ((37)), the objective function equation (36), and the initial estimation, the parameterized trajectory shape of the PE elements that satisfies the boundary conditions and thrust constraint is obtained by optimization.The complete process of PE-FFS can be expressed as follows: the PE elements are expressed by the Fourier series and discrete points are generated by the initial shape estimation equation (19) obtained from the polynomials. The initial estimate of all the coefficients are obtained from equations (17), (18), and ((20)–(22), and all the unknown coefficients are solved by the final optimization of the objective function in equation (36) and constraints in equations (28), (35), and (37).


## 2.1. Dynamic Modeling of a Two-Body System

A typical two-body problem [19] can be described as
(1)r¨+μr3r=f,where μ refers to the gravitational constant of the central body and f refers to the acceleration produced by the external force. In traditional shape-based methods, the state variable is represented by the polar radius r, the phase angle θ, and the height z along the z-axis. The components of equation (1) [32, 33] can be expressed as
(2)r¨−rθ̇+μrs3=fr,rθ¨+2ṙθ̇=fθ,z¨+μzs3=fz,where s=r2+z2 and fr,fθ,fz refer to the acceleration components in the local cylindrical coordinates. The definition of directions in the cylindrical coordinate system is shown in Figure 1.Figure 1
Diagram of the inertial and local cylindrical coordinate systems.The component in thez-axis periodically changes with the phase angle, and its peak and trough values are jointly determined by the semimajor axis, eccentricity, and inclination angle of the orbit. The schematic diagram of the component is shown in Figure 2.Figure 2
Schematic diagram of changes in thez-component during a multirevolution transfer.In the cylindrical coordinate system, if the inclination angle and eccentricity are small and there are few transfer cycles, an ordinary polynomial is suitable for the initial estimation of the component along thez-axis. As the inclination angle, eccentricity, and number of transfer cycles increase, convergence along the z-axis will become more difficult and the deviation in total velocity increment will increase gradually [31]. However, the orbit elements are slow variables to solve the problem. To avoid the singularity caused by zero eccentricity and zero inclination angle, the modified equinoctial elements are introduced and defined as follows [19]:
(3)p=a1−e2,f=ecosω+Ω,g=esinω+Ω,h=tani2cosΩ,k=tani2sinΩ,L=ω+Ω+θ,where a, e, i, ω, Ω, and θ denote the classical orbit elements. The two-body dynamic equation is expressed by the modified equinoctial elements in the Gaussian perturbation equation [19]:
(4)ṗ=2pp/μwf,ḟ=pμsinLfr+p/μw1+wcosL+fft−p/μwhsinL−kcosLgfn,ġ=−pμcosLfr+p/μw1+wsinL+gft+p/μwhsinL−kcosLffn,ḣ=s2pμcosL2wfn,k̇=s2pμsinL2wfn,L̇=p/μwhsinL−kcosLfn+μpw/p2,where s2=1+h2+k2, w=1+fcosL+gsinL, and fr,ft,fn denote the three acceleration components in the radial-transverse-normal coordinate system. The radial axis is the direction moving from the central gravitational body to the spacecraft, the normal axis is the direction of the trajectory angular momentum, and the transverse axis is determined by the right-hand rule. The components of the osculating orbit position rrx,ry,rz can be expressed by the modified equinoctial elements [19] as
(5)rx=rs2cosL+h2−k2cosL+2hksinL,ry=rs2sinL−h2−k2sinL+2hkcosL,rz=2rs2hsinL−kcosL,where r is the magnitude of the position vector, and r and s can be expressed as
(6)r=p1+fcosL+gsinL,s2=1+h2+k2.


## 2.2. Calculation of Guidance Law and On-Line Control Inputs

### 2.2.1. Guidance Law on Inverse Dynamics

One critical step of the shape-based design method is to represent the control variables with states. This process can be expressed in the form of an inverse model [34] when a typical system is considered:
(7)ẋ=Fx,t+Gx,tut,y=Hx,t.The derivative of the second term in equation (7) can be given by
(8)ẏ=∇Hx,tẋ=LfHx+LgHxut,where LfHx and LgHx, also called Lie derivatives [35], can be expressed as
(9)LfHx=∂Hx∂x⋅Fx,LgHx=∂Hx∂x⋅Gx.WhenLgHx≠0, it leads to the control variable by getting a single derivative:
(10)ut=LgHx−1ẏ−LfHx.Otherwise, it requiresk derivatives until LgLfk−1Hx≠0, so that
(11)ut=LgLfk−1Hx−1yk−LfkHx,where LfkHx and LgLfk−1Hx are defined, respectively, as
(12)Lfk−1Hx=LfLfk−2Hx=∂Lfk−2Hx∂xFx,LgLfk−1Hx=LgLfk−1Hx=∂Lfk−1Hx∂xGx.For the components of the cylindrical coordinatesr,θ,z, the first-order and second-order derivatives can be calculated; for the orbit elements, ṙ and r¨ are derived from equation (5). Derivatives of both dynamic models are then substituted into equation (11) to obtain the control variable ut.


## 2.2.1. Guidance Law on Inverse Dynamics

One critical step of the shape-based design method is to represent the control variables with states. This process can be expressed in the form of an inverse model [34] when a typical system is considered:
(7)ẋ=Fx,t+Gx,tut,y=Hx,t.The derivative of the second term in equation (7) can be given by
(8)ẏ=∇Hx,tẋ=LfHx+LgHxut,where LfHx and LgHx, also called Lie derivatives [35], can be expressed as
(9)LfHx=∂Hx∂x⋅Fx,LgHx=∂Hx∂x⋅Gx.WhenLgHx≠0, it leads to the control variable by getting a single derivative:
(10)ut=LgHx−1ẏ−LfHx.Otherwise, it requiresk derivatives until LgLfk−1Hx≠0, so that
(11)ut=LgLfk−1Hx−1yk−LfkHx,where LfkHx and LgLfk−1Hx are defined, respectively, as
(12)Lfk−1Hx=LfLfk−2Hx=∂Lfk−2Hx∂xFx,LgLfk−1Hx=LgLfk−1Hx=∂Lfk−1Hx∂xGx.For the components of the cylindrical coordinatesr,θ,z, the first-order and second-order derivatives can be calculated; for the orbit elements, ṙ and r¨ are derived from equation (5). Derivatives of both dynamic models are then substituted into equation (11) to obtain the control variable ut.


## 2.3. On-Line Control Based on Inverse Simulation

In a practical situation, the system is affected by external disturbances. A closed-loop control is required for the practical system. All the expected state variables are directly expressed with the independent time variable. Thus, the inverse simulation technique is quite suitable for the calculation of the practical control variables.The discrete system is(13)xtk+1=xtk+∫tktk+1ẋtk,uτdτ,ytk+1=Hxtk+1,where Hxtk+1 refers to the observation function. According to the current states obtained in the real-time system and the expected states, the expected control variable udtk can be obtained through iteration of the following [35]:
(14)un+1tk=untk−Jn−1Hxntk+1−ydtk+1,where J denotes the Jacobian matrix, which can be expressed as
(15)Jn=∂yerrorn∂u=yerrornun+δu−yerrornun−δu2δu,yerrornu=Hxntk+1−ydtk+1.The whole simulation framework is shown in Figure3.Figure 3
On-line control loop under the inverse simulation framework.The expected outputydtk+1 and the current state xtk are obtained based on the shaped trajectory and the measured value ytk, respectively. The expected control variable udtk is then calculated through the inverse simulation to form a closed loop.


## 2.4. Pseudoequinoctial Elements Based on the Finite Fourier Series

The state variable is represented in the form of a Fourier series:(16)λt=kλ2+∑n=1nλkaλncosnπTt+kbλnsinnπTt,where λ denotes the PE elements p~,f~,g~,h~,k~,L~ represented by the FFS, ka and kb represent the coefficients of each trigonometric term, nλ represents the number of the terms, and T denotes the TOF.


### 2.4.1. Calculation of the Initial Parameters

Based on the boundary conditions, the first two coefficients of the sine and cosine terms in each element expression can be calculated as [31]
(17)ka1=λ0−λf2−∑n=3nλkann≥3,n∈odd,ka2=λ0+λf−kλ2−∑n=3nλkann≥4,n∈even,(18)kb1=Tλ̇0−λ̇f2π−∑n=4nλnkbnn≥3,n∈odd,kb2=Tλ̇0+λ̇f4π−12∑n=4nλnkbnn≥4,n∈even,where λ0,λf denote the initial and final PE elements, respectively, and λ̇0,λ̇f represent the derivatives of the initial and final elements, respectively. For smooth transfer processes, the elements and their change rate at the beginning and end of the transfer orbit should be consistent with those of the initial and final orbits.After imposing the boundary conditions, the number of remaining unknown variables in equation (16) is 2nλ−3. The initial value of the unknown coefficients can be obtained from a series of collocation points that are in accordance with the initial estimation curves. The typical initial estimate can be obtained by a third-order polynomial that satisfies the boundary conditions:
(19)λinitt=λ̇0T+λ̇fT+2λ0−λftT3+3λf−λ0−2λ̇0T−λ̇fTtT2+λ̇0t+λ0.In this third-order polynomial curve, discrete points such as the Legendre-Gauss points and Legendre-Gauss-Lobatto points are selected as the collocation points.ndp refers to the number of points. By combining equations (16) and (19), the nonlinear equation set can be expressed as
(20)λinitndp×1=Trigndp×2nλ−3Kλ2nλ−3×1+Boundndp×1,where Kλ denotes the unknown Fourier series coefficients, and Trig and Bound are expressed as
(21)Trig=12cos3πTt⋯cosnλπTtsin3πTt⋯sinnλπTtnppr×1,Bound=ka1cosπTt+ka2cos2πTt+kb1sinπTt+kb2sin2πTtnppr×1.Eventually, the Moore-Penrose pseudoinverse [35] is applied to equation (20), through which the unknown coefficients Kλ can be expressed as
(22)Kλ=TrigTTrig−1TrigTλinit−Bound.


### 2.4.2. Practical Velocity and Acceleration Determination

After determining the initial values of all the coefficients, the next step is to calculate the thrust acceleration. According to the inverse simulation theory,ṙ and r¨ deduced from equation (5) can be transformed into the control acceleration u. The spacecraft velocity v=ṙ is expressed as
(23)v=∂r∂x∂x∂t+∂r∂L∂L∂t,where x refers to p,f,g,h,k. The osculating velocity vosc can be expressed as
(24)vosc=∂r∂LL̇osc=∂r∂Lμpr2.The actual velocity can also be expressed as(25)v=vosc+vgau,where vgau is the gauge velocity.By substituting equations (5) and (6) into equation (24), the osculating velocity component can be expressed as
(26)vx=−1p/μs2sinL+h2−k2sinL−2hkcosL+g−2fhk+h2−k2g,vy=−1p/μs2−cosL+h2−k2cosL+2hksinL−f+2ghk+h2−k2f,vz=2p/μs2hcosL+ksinL+fh+gk.When the dynamic equation satisfies the Gaussian perturbation equation, the osculating velocity and actual velocity coincide at each time, i.e., the osculating condition is satisfied andvgau≡0. The position and velocity can be obtained from equations (5) to (26), and vice versa. Thus, the modified equinoctial elements are equivalent to the osculating position and velocity. The PE elements, in the form of the Fourier series, are applied in order to ensure the bijection relationship between the motion states and pseudoelements. The position that satisfies equation (5) and the velocity at that time is
(27)v=∂r∂x~∂x~∂t+∂r∂L~∂L~∂t,where x~ refers to p~,f~,g~,h~,k~. In order to satisfy angular momentum conservation, an additional constraint is introduced and given by
(28)∂L~∂t=μp~r2.Thus, based on equations (5) and (27) and the hypothetical condition equation (28), the bijection relationship between the pseudoelements and contemporary motion states is obtained.In the calculation ofr¨, it is necessary to obtain a further time derivative of equation (24):
(29)r¨=∂vosc∂x~∂x~∂t+∂vosc∂L~∂L~∂t+∂vgau∂x~∂x~∂t+∂vgau∂L~∂L~∂t.For the modified equinoctial elements,vgau≡0, equation (29) will not provide any information on the external acceleration, i.e., r¨=−μ/r3r. However, it is impossible for pseudoelements to ensure that the gauge velocity is exactly zero. Equation (29) can be used to calculate the actual acceleration. The change in L~ with time is faster than in other elements. It can be deduced from equation (3) that when the external acceleration is lower than the local gravitational acceleration, the approximation condition is
(30)∂x~∂t≪∂L~∂t.As∂r/∂x~ and ∂r/∂L~ are of the same order, the pseudoelements should then satisfy
(31)∂r∂x~∂x~∂t≪∂r∂L~∂L~∂t.Then, the total velocity could be expressed asṙ≈∂r/∂L~∂L~/∂t and the osculating velocity and the position could be expressed as
(32)vosc=∂r∂LL̇osc≈∂r∂L~∂L~∂t,(33)r¨=∂v∂x~∂x~∂t+∂vgau∂L~L̇osc+∂vosc∂L~L̇osc.Thus, the osculating condition is basically satisfied. It can be inferred from equations (5), (27), and (33), according to the inverse simulation structure, that the thrust acceleration is expressed in the following form:
(34)F=Fλt,λ̇t,λ¨t.


### 2.4.3. Pseudoequinoctial Fourier Series Method

Eventually, within the transfer timet0−tf, the constraint and velocity increment can be separately expressed as
(35)c1F2−Fmax≤0,c2∂L∂t−L̇osc=0,(36)Δv=∫t0tfF2dt.The rectangular quadrature rule is utilized in this paper. In addition, to ensure a prograde trajectory, an additional constraint for the orbit direction can be added as(37)−c3L̇osc≤0,where c1−c3 denote the relaxation factors. The relationship between the fuel cost and final velocity increment is expressed as [18]
(38)mf=m01−exp−ΔvIspg0,where Isp denotes the specific impulse and g0 denotes the standard gravitational acceleration at sea level. According to the constraints of equations (35) and ((37)), the objective function equation (36), and the initial estimation, the parameterized trajectory shape of the PE elements that satisfies the boundary conditions and thrust constraint is obtained by optimization.The complete process of PE-FFS can be expressed as follows: the PE elements are expressed by the Fourier series and discrete points are generated by the initial shape estimation equation (19) obtained from the polynomials. The initial estimate of all the coefficients are obtained from equations (17), (18), and ((20)–(22), and all the unknown coefficients are solved by the final optimization of the objective function in equation (36) and constraints in equations (28), (35), and (37).


## 2.4.1. Calculation of the Initial Parameters

Based on the boundary conditions, the first two coefficients of the sine and cosine terms in each element expression can be calculated as [31]
(17)ka1=λ0−λf2−∑n=3nλkann≥3,n∈odd,ka2=λ0+λf−kλ2−∑n=3nλkann≥4,n∈even,(18)kb1=Tλ̇0−λ̇f2π−∑n=4nλnkbnn≥3,n∈odd,kb2=Tλ̇0+λ̇f4π−12∑n=4nλnkbnn≥4,n∈even,where λ0,λf denote the initial and final PE elements, respectively, and λ̇0,λ̇f represent the derivatives of the initial and final elements, respectively. For smooth transfer processes, the elements and their change rate at the beginning and end of the transfer orbit should be consistent with those of the initial and final orbits.After imposing the boundary conditions, the number of remaining unknown variables in equation (16) is 2nλ−3. The initial value of the unknown coefficients can be obtained from a series of collocation points that are in accordance with the initial estimation curves. The typical initial estimate can be obtained by a third-order polynomial that satisfies the boundary conditions:
(19)λinitt=λ̇0T+λ̇fT+2λ0−λftT3+3λf−λ0−2λ̇0T−λ̇fTtT2+λ̇0t+λ0.In this third-order polynomial curve, discrete points such as the Legendre-Gauss points and Legendre-Gauss-Lobatto points are selected as the collocation points.ndp refers to the number of points. By combining equations (16) and (19), the nonlinear equation set can be expressed as
(20)λinitndp×1=Trigndp×2nλ−3Kλ2nλ−3×1+Boundndp×1,where Kλ denotes the unknown Fourier series coefficients, and Trig and Bound are expressed as
(21)Trig=12cos3πTt⋯cosnλπTtsin3πTt⋯sinnλπTtnppr×1,Bound=ka1cosπTt+ka2cos2πTt+kb1sinπTt+kb2sin2πTtnppr×1.Eventually, the Moore-Penrose pseudoinverse [35] is applied to equation (20), through which the unknown coefficients Kλ can be expressed as
(22)Kλ=TrigTTrig−1TrigTλinit−Bound.


## 2.4.2. Practical Velocity and Acceleration Determination

After determining the initial values of all the coefficients, the next step is to calculate the thrust acceleration. According to the inverse simulation theory,ṙ and r¨ deduced from equation (5) can be transformed into the control acceleration u. The spacecraft velocity v=ṙ is expressed as
(23)v=∂r∂x∂x∂t+∂r∂L∂L∂t,where x refers to p,f,g,h,k. The osculating velocity vosc can be expressed as
(24)vosc=∂r∂LL̇osc=∂r∂Lμpr2.The actual velocity can also be expressed as(25)v=vosc+vgau,where vgau is the gauge velocity.By substituting equations (5) and (6) into equation (24), the osculating velocity component can be expressed as
(26)vx=−1p/μs2sinL+h2−k2sinL−2hkcosL+g−2fhk+h2−k2g,vy=−1p/μs2−cosL+h2−k2cosL+2hksinL−f+2ghk+h2−k2f,vz=2p/μs2hcosL+ksinL+fh+gk.When the dynamic equation satisfies the Gaussian perturbation equation, the osculating velocity and actual velocity coincide at each time, i.e., the osculating condition is satisfied andvgau≡0. The position and velocity can be obtained from equations (5) to (26), and vice versa. Thus, the modified equinoctial elements are equivalent to the osculating position and velocity. The PE elements, in the form of the Fourier series, are applied in order to ensure the bijection relationship between the motion states and pseudoelements. The position that satisfies equation (5) and the velocity at that time is
(27)v=∂r∂x~∂x~∂t+∂r∂L~∂L~∂t,where x~ refers to p~,f~,g~,h~,k~. In order to satisfy angular momentum conservation, an additional constraint is introduced and given by
(28)∂L~∂t=μp~r2.Thus, based on equations (5) and (27) and the hypothetical condition equation (28), the bijection relationship between the pseudoelements and contemporary motion states is obtained.In the calculation ofr¨, it is necessary to obtain a further time derivative of equation (24):
(29)r¨=∂vosc∂x~∂x~∂t+∂vosc∂L~∂L~∂t+∂vgau∂x~∂x~∂t+∂vgau∂L~∂L~∂t.For the modified equinoctial elements,vgau≡0, equation (29) will not provide any information on the external acceleration, i.e., r¨=−μ/r3r. However, it is impossible for pseudoelements to ensure that the gauge velocity is exactly zero. Equation (29) can be used to calculate the actual acceleration. The change in L~ with time is faster than in other elements. It can be deduced from equation (3) that when the external acceleration is lower than the local gravitational acceleration, the approximation condition is
(30)∂x~∂t≪∂L~∂t.As∂r/∂x~ and ∂r/∂L~ are of the same order, the pseudoelements should then satisfy
(31)∂r∂x~∂x~∂t≪∂r∂L~∂L~∂t.Then, the total velocity could be expressed asṙ≈∂r/∂L~∂L~/∂t and the osculating velocity and the position could be expressed as
(32)vosc=∂r∂LL̇osc≈∂r∂L~∂L~∂t,(33)r¨=∂v∂x~∂x~∂t+∂vgau∂L~L̇osc+∂vosc∂L~L̇osc.Thus, the osculating condition is basically satisfied. It can be inferred from equations (5), (27), and (33), according to the inverse simulation structure, that the thrust acceleration is expressed in the following form:
(34)F=Fλt,λ̇t,λ¨t.


## 2.4.3. Pseudoequinoctial Fourier Series Method

Eventually, within the transfer timet0−tf, the constraint and velocity increment can be separately expressed as
(35)c1F2−Fmax≤0,c2∂L∂t−L̇osc=0,(36)Δv=∫t0tfF2dt.The rectangular quadrature rule is utilized in this paper. In addition, to ensure a prograde trajectory, an additional constraint for the orbit direction can be added as(37)−c3L̇osc≤0,where c1−c3 denote the relaxation factors. The relationship between the fuel cost and final velocity increment is expressed as [18]
(38)mf=m01−exp−ΔvIspg0,where Isp denotes the specific impulse and g0 denotes the standard gravitational acceleration at sea level. According to the constraints of equations (35) and ((37)), the objective function equation (36), and the initial estimation, the parameterized trajectory shape of the PE elements that satisfies the boundary conditions and thrust constraint is obtained by optimization.The complete process of PE-FFS can be expressed as follows: the PE elements are expressed by the Fourier series and discrete points are generated by the initial shape estimation equation (19) obtained from the polynomials. The initial estimate of all the coefficients are obtained from equations (17), (18), and ((20)–(22), and all the unknown coefficients are solved by the final optimization of the objective function in equation (36) and constraints in equations (28), (35), and (37).


## 3. Two-Body Rendezvous Trajectory Design

In this section, the two-body rendezvous trajectory design using PE-FFS is investigated. It is assumed that in transfer problems including low Earth orbit (LEO), Earth-Mars, Earth-asteroid, and Earth-comet, the major flight is only affected by the central gravitational body. The differences in the order of magnitudes of each variable will make the computational numerical method less stable. Therefore, in this study, all the variables adopted are normalized. For the near-Earth space, the distance unit (DU) refers to the average radius of the Earth, and the2π time unit (TU) refers to the period of a circular orbit with the radius of the Earth. For the deep-space exploration problems, DU refers to the astronomical unit (AU), i.e., 1DU=149597870.66km. 2π TU equals the period of an orbit with the radius of 1 DU around the Sun, i.e., 1TU=58.135d. The pseudospectral method is completed by the general pseudospectral optimization software (GPOPS) [11]. The nonlinear optimization problem is solved by the sequential quadratic programming (SQP) algorithm in MATLAB’s nonlinear programming solver. A computer with a 3.4 GHz Intel Core i7-6700 processor, a Windows 7 operating system, and 8 GB RAM was used for the simulation.The data of the deep-space exploration missions reported in the literature [20] was used to study the trajectory design for three types of celestial bodies: planets, asteroids, and comets. The Kepler orbital elements of asteroid 1989ML and comet Tempel-1 and their epochs expressed in modified Julian date (MJD) are listed in Table 1. The orbital elements of Mars are calculated through the Jet Propulsion Laboratory (JPL) ephemeris [36]. In these missions, the maximum acceleration constraints were considered. The relaxation factors c1, c2, c3 were set as 10, 1, and 1, respectively, by trial and error; and the error tolerance was 1×10−2.Table 1
Classical orbital elements of 1989ML and Temple-1 at the given epoch MJD.
Parameter
a (DU)
e
i (°)
Ω (°)
ω (°)
M (°)
Epoch (MJD)
1989ML
1.2721
0.13649
4.3782
104.3571
183.3249
117.36689
53900
Tempel-1
3.14009
0.51159
10.5052
68.8818
179.3031
203.23760
56717During the calculation, the revolution number in 3D-FFS is traversed. Some initial determination was introduced to narrow the search range of the number of revolution. The relationship between the maximum accelerationFmax, phase angle θf, and transfer time Ttof, was considered. The phase angle can be expressed as
(39)θf=θ~+Nrev⋅2π,where θ~ represents the angle between the initial and final position vectors, which is typically defined as a positive value. If the angle is negative, then 2π needs to be added. Nrev refers to the number of orbit revolutions. If the instant semimajor axis of the transfer trajectory always lies in between those of the initial and final trajectories, then the maximum number of revolutions is
(40)Nrevmax=Ttof2πμa03−θ~2π.At the same time, considering the thruster capability, the minimum number of revolutions is(41)Nrevmin=maxTtof2πμaf3−θ~2π,v0−vf2πFmaxμa0+af/23−θ~2π,where v0 and vf denote the velocities of the initial and final trajectories, and a0 and af denote the semimajor axes of the initial and final trajectories, respectively. The physically impossible condition Nrevmax<Nrevmin indicates that the instant semimajor axis of the trajectory during the transfer process is slightly smaller than that of the initial trajectory or slightly larger than that of the final trajectory. In such a situation, Nrevmax and Nrevmin can be rounded instead of the upper and lower rounded-off numbers.


### 3.1. Earth to Mars

The search range in the literature [20], from January 1, 2020 to December 31, 2027, including four synodic periods, was adopted for the launch window analysis; the step size in the search space is 50 days. The TOF range changes from 500 to 2,000 days, and the step size in the search space is 50 days, with a maximum thrust acceleration Fmax of 1.5×10−4m/s2; the number of discrete points in each revolution is 15. All these points are used to enforce the constraints given by equations (35) and (37). The initial mass is 1,000 kg. The difference in the initial and final positions affects the trajectory shape, and the corresponding number of terms of PEs in equation (16) should also be changed. There are eight settings of the number of terms in the shape library, as presented in Table 2. During the calculations, only the optimal results in each of these eight settings were recorded.Table 2
Settings of the number of terms for the Earth-Mars rendezvous trajectory.
No.
1
2
3
4
5
6
7
8
np
4
4
3
4
3
3
3
3
nf=ng
3
3
3
2
4
4
3
2
nk=nh
2
2
3
2
2
2
2
2
nl
5
4
4
5
5
4
5
6The results are similar to those of previous studies [20, 22, 31]. Significant periodical patterns can be seen in Figure 4. The color bar refers to the velocity increment value, and the unit is in km/s. In general, more fuel is necessary for a long-term transfer and relatively reasonable trajectories are concentrated in the specific TOF within 1,200 days. The optimal solution yields a launch time of 9505 MJD2000, TOF of 950 days, a complete revolution of one, a velocity increment of 5.7684 km/s, and a corresponding shape number of five in Table 2. The trajectory and thrust acceleration are shown in Figures 5 and 6, respectively, and the comparison with the results of other methods is presented in Table 3.Figure 4
Launch window analysis in the Earth-Mars transfer mission.Figure 5
Earth-Mars rendezvous trajectory with the least amount of fuel.Figure 6
Trajectory acceleration profile with the least fuel requirement.Table 3
Comparison of the Earth-Mars rendezvous trajectory results with different methods.
No.
∆v (km/s)
F2max1×10−4 (m/s2)
Fmax (N)
PE-FFS
5.7684
1.226
0.15
3D-FFS [31]
5.7294
1.5
0.15
GPOPS [31]
5.7077
—
0.15
Indirect [31]
5.67
—
0.15
Spherical [20]
5.74
—
0.22
PE [20]
5.83
—
0.16
Hodographic [22]
5.77
1.5
—The above comparison indicates that the proposed method has the lowest acceleration peak of1.226×10−4m/s2. The velocity increment lies between those of PE and spherical methods and is better than that of the hodographic method. Although the different number of terms in Table 2 do not change much, it can guarantee that no infeasible points would appear in the chart. The results would be more obvious in the transfer cases with larger inclinations and eccentricities. The average computational time of PE-FFS and 3D-FFS for a single trajectory design is 0.64 and 0.7 s [31], respectively. The computational time of GPOPS is over 2 s.


### 3.2. Earth to Asteroid 1989ML

The synodic period of 1989ML with respect to the Earth is 3.3 years; thus, the same launch design range for Mars is employed with a step length of 50 days. If the TOF range changes from 400 to 1,000 days, then the step length is chosen as 10 days with the maximum thrust acceleration of3.1×10−4m/s2 and ten discrete points for each revolution. All these points should satisfy the constraints given by equations (35) and (37). The initial mass is 1,000 kg. The settings of the shape library are presented in Table 4. The launch window of the Earth to the asteroid 1989ML mission is shown in Figure 7. The meaning of the color bar is in accordance with that in Figure 4.Table 4
Setting of the terms for the rendezvous trajectory of asteroid 1989ML.
No.
1
2
3
4
5
6
np
4
3
4
3
4
3
nf=ng
3
2
2
4
2
4
nk=nh
2
2
2
2
2
3
nl
4
6
5
5
7
5Figure 7
Launch window analysis of the Earth to asteroid 1989ML mission.Because the number of terms is adjusted in the shape library, the changes in PE-FFS are smoother, and there are more feasible solutions in the whole process than those in 3D-FFS. The optimal results of the trajectory and thrust acceleration of the transfer are shown in Figures8 and 9, respectively.Figure 8
Optimal transfer trajectory of Earth-asteroid 1989ML.Figure 9
Profiles of the optimal thrust acceleration.The best launch window is 7955 MJD2000 with 440 days TOF and one complete revolution around the Sun. The corresponding shape number in Table4 is five, and the peak value of the thrust acceleration is 3.1×10−4m/s2. The comparison with other methods is presented in Table 5.Table 5
Comparison of the asteroid 1989ML rendezvous mission by different methods.
No.
∆v (km/s)
F2max1×10−4 (m/s2)
Fmax (N)
PE-FFS
4.58
3.1
0.30
3D-FFS [31]
4.20
3.1
0.30
GPOPS [31]
4.02
—
0.30
Indirect [31]
4.02
—
0.30
Spherical [20]
4.47
—
0.31
PE [20]
4.82
—
0.33
Hodographic [22]
4.22
4.7
—The average computational time of PE-FFS and 3D-FFS for a single trajectory design is 1.09 and 1.4 s [31], respectively. The computational time of GPOPS is over 4 s. The proposed method has a velocity increment of 4.58 km/s, and the results lie between those of the spherical method and the PE method. However, Figure 7 shows that broader feasible solutions can be obtained by this method compared to 3D-FFS in [31]. In particular, feasible solutions exist for every departure date with over 710 days TOF. Within the given launch time, feasible solutions can always be obtained with a velocity increment within 5.5 km/s, indicating that for the target rendezvous trajectory problem with a significant inclination angle and eccentricity (Table 1), the PE-FFS has a broader convergence field and provides more feasible solutions.


### 3.3. Earth to Comet Tempel-1

The range for the launch window is from January 1, 2000 to January 3, 2016, and the step length is 50 days. The TOF is from 400 to 1,500 days, and the step length is set as 30 days.Fmax is set as 7.1×10−4m/s2, and the initial mass is 1,000 kg. The range for the number of revolutions is 0–3, and the number of discrete points of each revolution is 20. All these points should satisfy the constraints given by equations (35) and (37). The setting of the terms in the shape library is presented in Table 6, and the results are shown in Figure 10. The color bar also refers to the velocity increment value in the unit of km/s.Table 6
Term setting of the Tempel-1 rendezvous trajectory.
No.
1
2
3
4
5
6
7
8
np
2
4
4
3
3
3
3
3
nf=ng
3
3
2
2
3
4
4
3
nk=nh
2
2
2
2
2
2
3
3
nl
5
5
6
4
4
4
4
5Figure 10
Launch window analysis of the comet Tempel-1 rendezvous mission.Figure10 shows that as transfer time decreases, there are fewer feasible solutions. For TOF over 800 days, the proposed method provides more feasible solutions than the spherical and PE methods [20]. Furthermore, within the entire design range, there are always feasible solutions for a velocity increment below 16.5 km/s. The curves of the optimal transfer trajectory and thrust acceleration are illustrated in Figures 11 and 12, respectively.Figure 11
Optimal transfer trajectory of Tempel-1 mission.Figure 12
Optimal Tempel-1 rendezvous trajectory acceleration.The best result for the entire process is 12.65 km/s for a departure epoch of 3550 MJD2000 with a TOF of 460 days and zero complete revolution around the Sun. The peak thrust acceleration is7.1e−4m/s2 and the corresponding shape number in Table 6 is eight. The comparison of results to those calculated using other methods is presented in Table 7.Table 7
Comparison of the comet Tempel-1 rendezvous trajectory results by different methods.
No.
Δv (km/s)
F2max1e−4 (m/s2)
Fmax (N)
PE-FFS
12.65
7.1
0.59
3D-FFS [31]
12.69
7.1
0.59
GPOPS [31]
11.45
—
0.59
Indirect [31]
11.45
—
0.59
Spherical [20]
11.13
—
1.4
PE [20]
13.44
—
1.13
Hodographic [22]
12.68
11.7
—The velocity increment of PE-FFS is better than that of 3D-FFS. At the same time, it is ensured that the thrust acceleration satisfies the constraint conditions. The average computational time of PE-FFS and 3D-FFS for a single trajectory design is 1.21 and 1.66 s [31]; both are less than that of GPOPS. This indicates that the proposed shape-based method is suitable for a large-design space research.From the analysis of the above near-Earth trajectory design and deep-space exploration trajectory design, it can be concluded that PE-FFS can be applied to design rendezvous trajectories from near-Earth to deep space, and the obtained thrust acceleration changes are smoother and slower than those obtained by 3D-FFS. The results of the optimal velocity increment are between those of the spherical and PE methods. Although the velocity increment is greater than that of 3D-FFS, with a gradual increase in the eccentricity and inclination angle of the transfer trajectory, more feasible solutions can be obtained than those obtained by 3D-FFS. For trajectory rendezvous problems with large eccentricity and inclination angle, improved velocity increment results can be obtained using PE-FFS. In addition, the introduction of a different number of terms in Tables2, 4, and 6 indicates that the numbers of the terms have changed slightly in different cases and the number of terms of each element is smaller than that in 3D-FFS. Although the fixed number of Fourier terms would gain results similar to those in [31], the change of the number could compensate for the infeasible area in the launch window chart. Thus, it is guaranteed that there are always feasible solutions that are acceptable for each search window.


## 3.1. Earth to Mars

The search range in the literature [20], from January 1, 2020 to December 31, 2027, including four synodic periods, was adopted for the launch window analysis; the step size in the search space is 50 days. The TOF range changes from 500 to 2,000 days, and the step size in the search space is 50 days, with a maximum thrust acceleration Fmax of 1.5×10−4m/s2; the number of discrete points in each revolution is 15. All these points are used to enforce the constraints given by equations (35) and (37). The initial mass is 1,000 kg. The difference in the initial and final positions affects the trajectory shape, and the corresponding number of terms of PEs in equation (16) should also be changed. There are eight settings of the number of terms in the shape library, as presented in Table 2. During the calculations, only the optimal results in each of these eight settings were recorded.Table 2
Settings of the number of terms for the Earth-Mars rendezvous trajectory.
No.
1
2
3
4
5
6
7
8
np
4
4
3
4
3
3
3
3
nf=ng
3
3
3
2
4
4
3
2
nk=nh
2
2
3
2
2
2
2
2
nl
5
4
4
5
5
4
5
6The results are similar to those of previous studies [20, 22, 31]. Significant periodical patterns can be seen in Figure 4. The color bar refers to the velocity increment value, and the unit is in km/s. In general, more fuel is necessary for a long-term transfer and relatively reasonable trajectories are concentrated in the specific TOF within 1,200 days. The optimal solution yields a launch time of 9505 MJD2000, TOF of 950 days, a complete revolution of one, a velocity increment of 5.7684 km/s, and a corresponding shape number of five in Table 2. The trajectory and thrust acceleration are shown in Figures 5 and 6, respectively, and the comparison with the results of other methods is presented in Table 3.Figure 4
Launch window analysis in the Earth-Mars transfer mission.Figure 5
Earth-Mars rendezvous trajectory with the least amount of fuel.Figure 6
Trajectory acceleration profile with the least fuel requirement.Table 3
Comparison of the Earth-Mars rendezvous trajectory results with different methods.
No.
∆v (km/s)
F2max1×10−4 (m/s2)
Fmax (N)
PE-FFS
5.7684
1.226
0.15
3D-FFS [31]
5.7294
1.5
0.15
GPOPS [31]
5.7077
—
0.15
Indirect [31]
5.67
—
0.15
Spherical [20]
5.74
—
0.22
PE [20]
5.83
—
0.16
Hodographic [22]
5.77
1.5
—The above comparison indicates that the proposed method has the lowest acceleration peak of1.226×10−4m/s2. The velocity increment lies between those of PE and spherical methods and is better than that of the hodographic method. Although the different number of terms in Table 2 do not change much, it can guarantee that no infeasible points would appear in the chart. The results would be more obvious in the transfer cases with larger inclinations and eccentricities. The average computational time of PE-FFS and 3D-FFS for a single trajectory design is 0.64 and 0.7 s [31], respectively. The computational time of GPOPS is over 2 s.


## 3.2. Earth to Asteroid 1989ML

The synodic period of 1989ML with respect to the Earth is 3.3 years; thus, the same launch design range for Mars is employed with a step length of 50 days. If the TOF range changes from 400 to 1,000 days, then the step length is chosen as 10 days with the maximum thrust acceleration of3.1×10−4m/s2 and ten discrete points for each revolution. All these points should satisfy the constraints given by equations (35) and (37). The initial mass is 1,000 kg. The settings of the shape library are presented in Table 4. The launch window of the Earth to the asteroid 1989ML mission is shown in Figure 7. The meaning of the color bar is in accordance with that in Figure 4.Table 4
Setting of the terms for the rendezvous trajectory of asteroid 1989ML.
No.
1
2
3
4
5
6
np
4
3
4
3
4
3
nf=ng
3
2
2
4
2
4
nk=nh
2
2
2
2
2
3
nl
4
6
5
5
7
5Figure 7
Launch window analysis of the Earth to asteroid 1989ML mission.Because the number of terms is adjusted in the shape library, the changes in PE-FFS are smoother, and there are more feasible solutions in the whole process than those in 3D-FFS. The optimal results of the trajectory and thrust acceleration of the transfer are shown in Figures8 and 9, respectively.Figure 8
Optimal transfer trajectory of Earth-asteroid 1989ML.Figure 9
Profiles of the optimal thrust acceleration.The best launch window is 7955 MJD2000 with 440 days TOF and one complete revolution around the Sun. The corresponding shape number in Table4 is five, and the peak value of the thrust acceleration is 3.1×10−4m/s2. The comparison with other methods is presented in Table 5.Table 5
Comparison of the asteroid 1989ML rendezvous mission by different methods.
No.
∆v (km/s)
F2max1×10−4 (m/s2)
Fmax (N)
PE-FFS
4.58
3.1
0.30
3D-FFS [31]
4.20
3.1
0.30
GPOPS [31]
4.02
—
0.30
Indirect [31]
4.02
—
0.30
Spherical [20]
4.47
—
0.31
PE [20]
4.82
—
0.33
Hodographic [22]
4.22
4.7
—The average computational time of PE-FFS and 3D-FFS for a single trajectory design is 1.09 and 1.4 s [31], respectively. The computational time of GPOPS is over 4 s. The proposed method has a velocity increment of 4.58 km/s, and the results lie between those of the spherical method and the PE method. However, Figure 7 shows that broader feasible solutions can be obtained by this method compared to 3D-FFS in [31]. In particular, feasible solutions exist for every departure date with over 710 days TOF. Within the given launch time, feasible solutions can always be obtained with a velocity increment within 5.5 km/s, indicating that for the target rendezvous trajectory problem with a significant inclination angle and eccentricity (Table 1), the PE-FFS has a broader convergence field and provides more feasible solutions.


## 3.3. Earth to Comet Tempel-1

The range for the launch window is from January 1, 2000 to January 3, 2016, and the step length is 50 days. The TOF is from 400 to 1,500 days, and the step length is set as 30 days.Fmax is set as 7.1×10−4m/s2, and the initial mass is 1,000 kg. The range for the number of revolutions is 0–3, and the number of discrete points of each revolution is 20. All these points should satisfy the constraints given by equations (35) and (37). The setting of the terms in the shape library is presented in Table 6, and the results are shown in Figure 10. The color bar also refers to the velocity increment value in the unit of km/s.Table 6
Term setting of the Tempel-1 rendezvous trajectory.
No.
1
2
3
4
5
6
7
8
np
2
4
4
3
3
3
3
3
nf=ng
3
3
2
2
3
4
4
3
nk=nh
2
2
2
2
2
2
3
3
nl
5
5
6
4
4
4
4
5Figure 10
Launch window analysis of the comet Tempel-1 rendezvous mission.Figure10 shows that as transfer time decreases, there are fewer feasible solutions. For TOF over 800 days, the proposed method provides more feasible solutions than the spherical and PE methods [20]. Furthermore, within the entire design range, there are always feasible solutions for a velocity increment below 16.5 km/s. The curves of the optimal transfer trajectory and thrust acceleration are illustrated in Figures 11 and 12, respectively.Figure 11
Optimal transfer trajectory of Tempel-1 mission.Figure 12
Optimal Tempel-1 rendezvous trajectory acceleration.The best result for the entire process is 12.65 km/s for a departure epoch of 3550 MJD2000 with a TOF of 460 days and zero complete revolution around the Sun. The peak thrust acceleration is7.1e−4m/s2 and the corresponding shape number in Table 6 is eight. The comparison of results to those calculated using other methods is presented in Table 7.Table 7
Comparison of the comet Tempel-1 rendezvous trajectory results by different methods.
No.
Δv (km/s)
F2max1e−4 (m/s2)
Fmax (N)
PE-FFS
12.65
7.1
0.59
3D-FFS [31]
12.69
7.1
0.59
GPOPS [31]
11.45
—
0.59
Indirect [31]
11.45
—
0.59
Spherical [20]
11.13
—
1.4
PE [20]
13.44
—
1.13
Hodographic [22]
12.68
11.7
—The velocity increment of PE-FFS is better than that of 3D-FFS. At the same time, it is ensured that the thrust acceleration satisfies the constraint conditions. The average computational time of PE-FFS and 3D-FFS for a single trajectory design is 1.21 and 1.66 s [31]; both are less than that of GPOPS. This indicates that the proposed shape-based method is suitable for a large-design space research.From the analysis of the above near-Earth trajectory design and deep-space exploration trajectory design, it can be concluded that PE-FFS can be applied to design rendezvous trajectories from near-Earth to deep space, and the obtained thrust acceleration changes are smoother and slower than those obtained by 3D-FFS. The results of the optimal velocity increment are between those of the spherical and PE methods. Although the velocity increment is greater than that of 3D-FFS, with a gradual increase in the eccentricity and inclination angle of the transfer trajectory, more feasible solutions can be obtained than those obtained by 3D-FFS. For trajectory rendezvous problems with large eccentricity and inclination angle, improved velocity increment results can be obtained using PE-FFS. In addition, the introduction of a different number of terms in Tables2, 4, and 6 indicates that the numbers of the terms have changed slightly in different cases and the number of terms of each element is smaller than that in 3D-FFS. Although the fixed number of Fourier terms would gain results similar to those in [31], the change of the number could compensate for the infeasible area in the launch window chart. Thus, it is guaranteed that there are always feasible solutions that are acceptable for each search window.


## 4. Low-Thrust Trajectory Design and Analysis for CRTBP

In the above cases, transfer processes are frequently affected by the central gravitational body; hence, the two-body assumption was applied. When the period under the influence of a third gravitational body cannot be neglected, it is necessary to consider a new method for trajectory design. Using the Earth-Moon rendezvous trajectory design as an example, the CRTBP, a method that employs two two-body models to calculate the three-body low-thrust trajectory, is proposed. The free intermediate section without low-thrust acceleration during the Earth-Moon transfer was omitted from this study to shorten the transfer time. The capture energy in the LSOI was first analyzed, and the results were used as the final criterion for the lunar capture. The geocentric trajectory of the Earth-Moon transfer was designed using PE-FFS and compared to the 3D-FFS method.


### 4.1. Analysis of the Capture Energy in Lunar Sphere of Influence

According to Laplace’s definition of LSOI, the radius of the LSOI is [37]
(42)rm=re−mmmme2.5,where re−m denotes the Earth-Moon distance, and mm and me denote the mass of the Moon and Earth, respectively. The motion of the vehicle along the trajectory entering the LSOI was analyzed, as shown in Figure 13. The capture angle εL is the angle between the velocity vector vBL and the position vector rBL at capture point B, where rBL=rm. The position and velocity vectors at perilune C are denoted as rCL and vCL, respectively, and the phase angle is θ. Thus, the motion in the orbit plane can be represented by the polar coordinate system satisfying the first two terms of equation (2), where s degenerates to r.Figure 13
Diagram of the selenocentric trajectory.When a Kepler orbit is considered, the semilatus rectum of the trajectorypL can be calculated using the critical capture parameters vBL and εL. The semimajor axis aL, eccentricity eL, and phase angle θ are expressed as [37]
(43)aL=rm2−γBL,pL=rmγBLsin2εL,eL=1−pLaL,θ=arccos1eLpLrm−1,where γBL denotes the energy factor, which can be expressed as
(44)γBL=rmvBL2μL,where μL denotes the lunar gravitational constant. The perilune radius rCL can be obtained from equation (43). When rCL is less than the radius of the Moon, the vehicle will hit the Moon. For the hyperbolic and parabolic trajectories [37], the transfer time from point B to point C is
(45)ΔtL=−aL3μLeLsinhHBL−HBL,ΔtL=aL3μLEBL−eLsinEBL,where
(46)HBL=2arctaneL−1eL+1tanθ2,EBL=2arctan1−eLeL+1tanθ2.Except for the case of a hard landing, two capture criteria are given:Criterion 1.
During the transfer process, the range of feasible radiusr is between the lunar radius and the boundary of the LSOI.Criterion 2.
During the transfer process, the velocity is lower than the local escape velocity, which is expressed as(47)v<vescape=2μLr.When the initial state of the vehicle is propagated without external thrust, the relationship between the capture parametersvBL and εL and perilune distance rCL, and that between vBL and εL and transfer time ΔtL can be displayed in Figures 14 and 15, respectively. The color bars refer to the times about perilune distance to the lunar radius and transfer time in seconds, respectively. The perilune distance was normalized using the lunar radius as the unit distance. The vacancy area of the figures does not satisfy Criterion 1.Figure 14
rCL as a function of vBL and εL.Figure 15
ΔtL as a function of vBL and εL.The above figures indicate that ifεL is within the range of -0.1 to 0.1 or vBL is lower than 70 m/s, then the vehicle impacts the Moon. When the velocity is higher than 385 m/s, the escape velocity at the gravitational boundary, and εL is larger, the vehicle will fly out of the LSOI during the transfer time ΔtL. As vBL increases, the feasible range of εL shrinks gradually. Within the range of vBL larger than 385 m/s, the rCL and transfer time can be feasible only under the condition of an extremely small εL, but the vehicle does not satisfy Criterion 2 and will fly out of the LSOI after 2ΔtL.To expand the area satisfying Criterion2, thrust is applied to brake the vehicle. Constant thrust acceleration from 1×10−6m/s2 to 1×10−1m/s2 is imposed opposite the direction of the velocity within a time of 2ΔtL. The relationship between the minimum acceleration to satisfy two capture criteria and the capture parameters is illustrated Figure 16, in which the contour lines refer to the exponents of the acceleration. The relationship between the perilune distance under the deceleration and capture parameters is depicted in Figure 17. As before, the color bar refers to the times of the lunar radius.Figure 16
Fmax as a function of vBL and εL.Figure 17
rCL as a function of vBL and εL.From Figures16 and 17, it can be inferred that the feasible range enlarges due to the deceleration thrust. The dark-colored area indicates that the vehicle is automatically captured. The light-colored area indicates that the vehicle is captured by the lunar gravity under the thruster brake. The area between -3 and -2 indicates that with increasing vBL, the thrust for the lunar capture has to increase. rCL decreases with smaller εL. In the area between -2 and -1, rCL is at the edge of the LSOI and the amplitude of deceleration is large.


### 4.2. Earth-Moon Transfer Launch Window Design

The geocentric two-body model was adopted in the escape section, and the proposed PE-FFS was applied. With the Moon being the rendezvous target, the minimum number of transfer revolutions is determined by the revolution formula of the tangential thrust escape trajectory [13], which can be expressed as
(48)Nrevmin=g08πamax.During the transfer process, the trajectory will pass through the LSOI. At that moment, the positionrBL and velocity vBL relative to the Moon can be expressed as
(49)rBL=r−rL,vBL=v−vL.The angleεL between the position and the velocity can be expressed as
(50)εL=arcsinrBL×vBLrBLvBL.The range of the angle is between−π/2 and π/2.The first step is to design the whole trajectory with the initial GEO states and the final motion states of the Moon by PE-FFS under the two-body dynamics. Because the tree-body dynamics is considered, only the segment out of the LSOI can be used as the geocentric trajectory. The moment of reaching the LSOI can be calculated by performing an iteration with equation (5). The capture parameters can then be calculated using equations (26), (49), and (50) to further determine the feasibility of the trajectory.The parameters of the first case in [29] were employed for the simulation analysis. The vehicle transfers from the geosynchronous orbit (GEO) to the lunar trajectory. The maximum acceleration is 3.237×10−2m/s2 and the range of transfer time is 5 to 8.5 days. It starts at different moments within the GEO trajectory period. The number of terms is set as np=2or3,nf=ng=2,nk=nh=2,andnl=6. The number of collocation points per revolution is 4.2. It can be derived from equation (48) that the minimum number of revolutions is 12 and the traversal range is from 12 to 22. The capture parameters are calculated according to equations (49) and (50), which are presented in Table 8. Δvout refers to the velocity increment from the GEO to the LSOI. Δv refers to the total velocity increment under the two-body dynamics. The proportion is the ratio of the Δvout to the Δv.Table 8
Capture parameters of the Earth-Moon transfer trajectory.
Transfer time (days)
Capture time (days)
Rev no.
εL (rad)
vBL (m/s)
Δvout (km/s)
Δv (km/s)
Proportion (%)
5
4.336
13
0.874
170.5
7.483
8.552
87.50
5.5
4.779
14
0.930
173.5
7.366
8.299
88.76
6
5.271
15
1.049
175.5
7.306
8.087
90.34
6.5
5.779
16
1.142
177.1
7.265
7.932
91.59
7
6.256
18
1.205
180.1
7.275
7.857
92.59
7.5
6.760
19
1.265
180.0
7.203
7.702
93.52
8
7.267
20
1.331
180.9
7.215
7.628
94.59
8.5
7.766
21
1.356
181.0
7.131
7.503
95.04Within the transfer time of 5–8.5 days, the initial phase of the GEO trajectory has slight effects on the velocity increment, and the velocity increment decreases as the TOF increases. The data in the above table shows that the ratio of the velocity increment outside the LSOI to the total velocity increases when the transfer time is prolonged and all the capture angles are within 0.5–1.57 rad. The trajectory appears to be retrograding around the Moon in the velocity range of 170–180 m/s. The results from literature [29, 38] are listed in Table 9. The velocity increments are the results for the entire process, including the escape phase, intermediate phase, and capture phase, whereas the cycle and transfer times are the results for the escape and intermediate phases. The trajectory direction means the orbit direction of the final phase around the Moon. The former two phases are outside the LSOI, while the last phase is in the LSOI.Table 9
Comparison of results of different methods [29, 38].
Method
Maximum thrust (1×10−2m/s2)
Transfer cycle
Transfer time (days)
Trajectory direction
Δv (km/s)
CPU (s)
Three-stage approach [38]
2.943
12
6.81
Prograde
7.05
—
2D-FFS [29]
3.237
12
6.45
Prograde
8.14
129The velocity increment for the escape and intermediate phases in [29] is approximately 7 km/s. The closest transfer time is 6.76 days in Table 8. The velocity increment of the escape and intermediate phases is 7.203 km/s, which is 2.9% larger than that of 2D-FFS, due to the deficiency of the intermediate phase. However, the computational time of PE-FFS is 7.35 s, which is suitable for the initial estimation and launch window design. Comparing capture parameters with the results in Figures 14–17, the low-thrust deceleration ensures that the vehicle is captured by the Moon; hence, the result can be regarded as the initial solution for the geocentric transfer trajectory. The trajectory and thrust acceleration curves obtained from PE-FFS are shown in Figures 18 and 19, respectively.Figure 18
Geocentric section of Earth-Moon transfer trajectory.Figure 19
Profile of the thrust acceleration.Within 1.3 days after the start of the transfer, the vehicle accelerates along the Earth’s orbit and the peak acceleration reaches 0.03 m/s2. Subsequently, as the acceleration decreases, the vehicle flies toward the Moon. After 3.5 days, the thrust descends below 30% of the acceleration peak. This phase can be regarded as the intermediate phase. After obtaining the geocentric section in the Earth-Moon transfer, the position and velocity of the vehicle relative to the Moon at the capture point can be used to further design the selenocentric trajectory. The Moon becomes the central gravitational body and the PE-FFS can still be applied.


## 4.1. Analysis of the Capture Energy in Lunar Sphere of Influence

According to Laplace’s definition of LSOI, the radius of the LSOI is [37]
(42)rm=re−mmmme2.5,where re−m denotes the Earth-Moon distance, and mm and me denote the mass of the Moon and Earth, respectively. The motion of the vehicle along the trajectory entering the LSOI was analyzed, as shown in Figure 13. The capture angle εL is the angle between the velocity vector vBL and the position vector rBL at capture point B, where rBL=rm. The position and velocity vectors at perilune C are denoted as rCL and vCL, respectively, and the phase angle is θ. Thus, the motion in the orbit plane can be represented by the polar coordinate system satisfying the first two terms of equation (2), where s degenerates to r.Figure 13
Diagram of the selenocentric trajectory.When a Kepler orbit is considered, the semilatus rectum of the trajectorypL can be calculated using the critical capture parameters vBL and εL. The semimajor axis aL, eccentricity eL, and phase angle θ are expressed as [37]
(43)aL=rm2−γBL,pL=rmγBLsin2εL,eL=1−pLaL,θ=arccos1eLpLrm−1,where γBL denotes the energy factor, which can be expressed as
(44)γBL=rmvBL2μL,where μL denotes the lunar gravitational constant. The perilune radius rCL can be obtained from equation (43). When rCL is less than the radius of the Moon, the vehicle will hit the Moon. For the hyperbolic and parabolic trajectories [37], the transfer time from point B to point C is
(45)ΔtL=−aL3μLeLsinhHBL−HBL,ΔtL=aL3μLEBL−eLsinEBL,where
(46)HBL=2arctaneL−1eL+1tanθ2,EBL=2arctan1−eLeL+1tanθ2.Except for the case of a hard landing, two capture criteria are given:Criterion 1.
During the transfer process, the range of feasible radiusr is between the lunar radius and the boundary of the LSOI.Criterion 2.
During the transfer process, the velocity is lower than the local escape velocity, which is expressed as(47)v<vescape=2μLr.When the initial state of the vehicle is propagated without external thrust, the relationship between the capture parametersvBL and εL and perilune distance rCL, and that between vBL and εL and transfer time ΔtL can be displayed in Figures 14 and 15, respectively. The color bars refer to the times about perilune distance to the lunar radius and transfer time in seconds, respectively. The perilune distance was normalized using the lunar radius as the unit distance. The vacancy area of the figures does not satisfy Criterion 1.Figure 14
rCL as a function of vBL and εL.Figure 15
ΔtL as a function of vBL and εL.The above figures indicate that ifεL is within the range of -0.1 to 0.1 or vBL is lower than 70 m/s, then the vehicle impacts the Moon. When the velocity is higher than 385 m/s, the escape velocity at the gravitational boundary, and εL is larger, the vehicle will fly out of the LSOI during the transfer time ΔtL. As vBL increases, the feasible range of εL shrinks gradually. Within the range of vBL larger than 385 m/s, the rCL and transfer time can be feasible only under the condition of an extremely small εL, but the vehicle does not satisfy Criterion 2 and will fly out of the LSOI after 2ΔtL.To expand the area satisfying Criterion2, thrust is applied to brake the vehicle. Constant thrust acceleration from 1×10−6m/s2 to 1×10−1m/s2 is imposed opposite the direction of the velocity within a time of 2ΔtL. The relationship between the minimum acceleration to satisfy two capture criteria and the capture parameters is illustrated Figure 16, in which the contour lines refer to the exponents of the acceleration. The relationship between the perilune distance under the deceleration and capture parameters is depicted in Figure 17. As before, the color bar refers to the times of the lunar radius.Figure 16
Fmax as a function of vBL and εL.Figure 17
rCL as a function of vBL and εL.From Figures16 and 17, it can be inferred that the feasible range enlarges due to the deceleration thrust. The dark-colored area indicates that the vehicle is automatically captured. The light-colored area indicates that the vehicle is captured by the lunar gravity under the thruster brake. The area between -3 and -2 indicates that with increasing vBL, the thrust for the lunar capture has to increase. rCL decreases with smaller εL. In the area between -2 and -1, rCL is at the edge of the LSOI and the amplitude of deceleration is large.


## 4.2. Earth-Moon Transfer Launch Window Design

The geocentric two-body model was adopted in the escape section, and the proposed PE-FFS was applied. With the Moon being the rendezvous target, the minimum number of transfer revolutions is determined by the revolution formula of the tangential thrust escape trajectory [13], which can be expressed as
(48)Nrevmin=g08πamax.During the transfer process, the trajectory will pass through the LSOI. At that moment, the positionrBL and velocity vBL relative to the Moon can be expressed as
(49)rBL=r−rL,vBL=v−vL.The angleεL between the position and the velocity can be expressed as
(50)εL=arcsinrBL×vBLrBLvBL.The range of the angle is between−π/2 and π/2.The first step is to design the whole trajectory with the initial GEO states and the final motion states of the Moon by PE-FFS under the two-body dynamics. Because the tree-body dynamics is considered, only the segment out of the LSOI can be used as the geocentric trajectory. The moment of reaching the LSOI can be calculated by performing an iteration with equation (5). The capture parameters can then be calculated using equations (26), (49), and (50) to further determine the feasibility of the trajectory.The parameters of the first case in [29] were employed for the simulation analysis. The vehicle transfers from the geosynchronous orbit (GEO) to the lunar trajectory. The maximum acceleration is 3.237×10−2m/s2 and the range of transfer time is 5 to 8.5 days. It starts at different moments within the GEO trajectory period. The number of terms is set as np=2or3,nf=ng=2,nk=nh=2,andnl=6. The number of collocation points per revolution is 4.2. It can be derived from equation (48) that the minimum number of revolutions is 12 and the traversal range is from 12 to 22. The capture parameters are calculated according to equations (49) and (50), which are presented in Table 8. Δvout refers to the velocity increment from the GEO to the LSOI. Δv refers to the total velocity increment under the two-body dynamics. The proportion is the ratio of the Δvout to the Δv.Table 8
Capture parameters of the Earth-Moon transfer trajectory.
Transfer time (days)
Capture time (days)
Rev no.
εL (rad)
vBL (m/s)
Δvout (km/s)
Δv (km/s)
Proportion (%)
5
4.336
13
0.874
170.5
7.483
8.552
87.50
5.5
4.779
14
0.930
173.5
7.366
8.299
88.76
6
5.271
15
1.049
175.5
7.306
8.087
90.34
6.5
5.779
16
1.142
177.1
7.265
7.932
91.59
7
6.256
18
1.205
180.1
7.275
7.857
92.59
7.5
6.760
19
1.265
180.0
7.203
7.702
93.52
8
7.267
20
1.331
180.9
7.215
7.628
94.59
8.5
7.766
21
1.356
181.0
7.131
7.503
95.04Within the transfer time of 5–8.5 days, the initial phase of the GEO trajectory has slight effects on the velocity increment, and the velocity increment decreases as the TOF increases. The data in the above table shows that the ratio of the velocity increment outside the LSOI to the total velocity increases when the transfer time is prolonged and all the capture angles are within 0.5–1.57 rad. The trajectory appears to be retrograding around the Moon in the velocity range of 170–180 m/s. The results from literature [29, 38] are listed in Table 9. The velocity increments are the results for the entire process, including the escape phase, intermediate phase, and capture phase, whereas the cycle and transfer times are the results for the escape and intermediate phases. The trajectory direction means the orbit direction of the final phase around the Moon. The former two phases are outside the LSOI, while the last phase is in the LSOI.Table 9
Comparison of results of different methods [29, 38].
Method
Maximum thrust (1×10−2m/s2)
Transfer cycle
Transfer time (days)
Trajectory direction
Δv (km/s)
CPU (s)
Three-stage approach [38]
2.943
12
6.81
Prograde
7.05
—
2D-FFS [29]
3.237
12
6.45
Prograde
8.14
129The velocity increment for the escape and intermediate phases in [29] is approximately 7 km/s. The closest transfer time is 6.76 days in Table 8. The velocity increment of the escape and intermediate phases is 7.203 km/s, which is 2.9% larger than that of 2D-FFS, due to the deficiency of the intermediate phase. However, the computational time of PE-FFS is 7.35 s, which is suitable for the initial estimation and launch window design. Comparing capture parameters with the results in Figures 14–17, the low-thrust deceleration ensures that the vehicle is captured by the Moon; hence, the result can be regarded as the initial solution for the geocentric transfer trajectory. The trajectory and thrust acceleration curves obtained from PE-FFS are shown in Figures 18 and 19, respectively.Figure 18
Geocentric section of Earth-Moon transfer trajectory.Figure 19
Profile of the thrust acceleration.Within 1.3 days after the start of the transfer, the vehicle accelerates along the Earth’s orbit and the peak acceleration reaches 0.03 m/s2. Subsequently, as the acceleration decreases, the vehicle flies toward the Moon. After 3.5 days, the thrust descends below 30% of the acceleration peak. This phase can be regarded as the intermediate phase. After obtaining the geocentric section in the Earth-Moon transfer, the position and velocity of the vehicle relative to the Moon at the capture point can be used to further design the selenocentric trajectory. The Moon becomes the central gravitational body and the PE-FFS can still be applied.


## 5. Conclusions

In this study, a new method for the shape-based trajectory design was proposed. The inverse dynamic model for the two-body system was established and the relationship between acceleration and the pseudoequinoctial elements was proposed. Following this, the proposed method was examined using two-body system problems for deep-space explorations such as Earth-Mars, Earth-1989ML, and Earth-Tempel-1. Under the conditions of small inclination angle and eccentricity, the velocity increments of the proposed method are between those of the spherical parameter method and pseudoelement method. With a large inclination angle and eccentricity, the proposed method yields more feasible solutions than the finite Fourier series method. For three-body trajectory transfer problems, the proposed method combined with the patched conics method was presented to analyze the low-thrust Earth-Moon transfer window. The combined method was applied to comprehensively design the initial low-thrust geocentric trajectory for Earth-Moon transfer, which can improve the design efficiency and ensure precision at the same time.



---

*Source: 1364834-2019-11-04.xml*