Research Article  Open Access
Wanmeng Zhou, Haiyang Li, Hua Wang, Ran Ding, "LowThrust Trajectory Design Using Finite Fourier Series Approximation of Pseudoequinoctial Elements", International Journal of Aerospace Engineering, vol. 2019, Article ID 1364834, 18 pages, 2019. https://doi.org/10.1155/2019/1364834
LowThrust Trajectory Design Using Finite Fourier Series Approximation of Pseudoequinoctial Elements
Abstract
A new lowthrust 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 online 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 twobody problems, three cases were studied: the Earth to Mars, 1989ML, and Tempel1 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 threebody problem was solved for the first time using the pseudoequinoctial finite Fourier series method combined with the patched conics method. The lowthrust EarthMoon transfer was analyzed, and the results show that this method improves window analysis efficiency and guarantees precision of the initial geocentric trajectory for the lowthrust transfer.
1. Introduction
In recent years, lowthrust spacecraft have increasingly attracted interest in the field of engineering because of their longlasting and fuelsaving features. These include the National Aeronautics and Space Administration’s (NASA) Deep Space 1 [1] and Dawn [2], the European Space Agency’s (ESA) SMART1 [3], and BepiColombo [4], a joint mission between ESA and Japan Aerospace Exploration Agency (JAXA). The trajectory design of these lowthrust vehicles has also attracted much attention. In a global trajectory optimization competition (GTOC 2), the trajectory design problem for lowthrust vehicles was proposed to explore asteroids and issues associated with them [5]. The challenges in the search for the optimal trajectory of a lowthrust, deepspace exploration vehicle are greater than those of an impulsivethrust 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 twopoint 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 shapebased 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 shapebased methods on twodimensional 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 threedimensional (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 switchingtype function with sufficient trigonometric function terms [28], which introduces more shape categories into the trajectory design. On the other hand, most research on the shapebased methods is applied on twobody problems, but few are applied for circular restricted threebody problems (CRTBPs). FFS combined with the simplified analysis method is the only shapebased 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 nearoptimum 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 3DFFS in a cylindrical coordinate system has limited approximation capability. During the longtime transfer, the motion in the 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 PEFFS. The PEs are similar to modified equinoctial elements (MEEs), which can provide smooth variations of instantaneous motion states compared to the axis components of a cylindrical coordinate system [32]. For twobody 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 PEFFS 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 twobody 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 PEFFS method. The third section presents the results of studies of the twobody rendezvous cases, which include transfers to Mars, an asteroid, and a comet. The results from the PEFFS method are compared with those from the other methods. Finally, the CRTBP for the EarthMoon transfer is investigated. The quick trajectory design and the launch window analysis for the EarthMoon mission are realized by combining the PEFFS and patched conics method.
2. Pseudoequinoctial Finite Fourier Series Method
2.1. Dynamic Modeling of a TwoBody System
A typical twobody problem [19] can be described as where refers to the gravitational constant of the central body and refers to the acceleration produced by the external force. In traditional shapebased methods, the state variable is represented by the polar radius , the phase angle , and the height along the axis. The components of equation (1) [32, 33] can be expressed as where and refer to the acceleration components in the local cylindrical coordinates. The definition of directions in the cylindrical coordinate system is shown in Figure 1.
The component in the 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.
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 the axis. As the inclination angle, eccentricity, and number of transfer cycles increase, convergence along the 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]: where , , , , , and denote the classical orbit elements. The twobody dynamic equation is expressed by the modified equinoctial elements in the Gaussian perturbation equation [19]: where , , and denote the three acceleration components in the radialtransversenormal 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 righthand rule. The components of the osculating orbit position can be expressed by the modified equinoctial elements [19] as where is the magnitude of the position vector, and and can be expressed as
2.2. Calculation of Guidance Law and OnLine Control Inputs
2.2.1. Guidance Law on Inverse Dynamics
One critical step of the shapebased 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:
The derivative of the second term in equation (7) can be given by where and , also called Lie derivatives [35], can be expressed as
When , it leads to the control variable by getting a single derivative:
Otherwise, it requires derivatives until , so that where and are defined, respectively, as
For the components of the cylindrical coordinates , the firstorder and secondorder derivatives can be calculated; for the orbit elements, and are derived from equation (5). Derivatives of both dynamic models are then substituted into equation (11) to obtain the control variable .
2.3. OnLine Control Based on Inverse Simulation
In a practical situation, the system is affected by external disturbances. A closedloop 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 where refers to the observation function. According to the current states obtained in the realtime system and the expected states, the expected control variable can be obtained through iteration of the following [35]: where denotes the Jacobian matrix, which can be expressed as
The whole simulation framework is shown in Figure 3.
The expected output and the current state are obtained based on the shaped trajectory and the measured value , respectively. The expected control variable 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: where denotes the PE elements represented by the FFS, and represent the coefficients of each trigonometric term, represents the number of the terms, and 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] where denote the initial and final PE elements, respectively, and 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 . 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 thirdorder polynomial that satisfies the boundary conditions:
In this thirdorder polynomial curve, discrete points such as the LegendreGauss points and LegendreGaussLobatto points are selected as the collocation points. refers to the number of points. By combining equations (16) and (19), the nonlinear equation set can be expressed as where denotes the unknown Fourier series coefficients, and and are expressed as
Eventually, the MoorePenrose pseudoinverse [35] is applied to equation (20), through which the unknown coefficients can be expressed as
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, and deduced from equation (5) can be transformed into the control acceleration . The spacecraft velocity is expressed as where refers to . The osculating velocity can be expressed as
The actual velocity can also be expressed as where is the gauge velocity.
By substituting equations (5) and (6) into equation (24), the osculating velocity component can be expressed as
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 and . 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 where refers to . In order to satisfy angular momentum conservation, an additional constraint is introduced and given by
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 of , it is necessary to obtain a further time derivative of equation (24):
For the modified equinoctial elements, , equation (29) will not provide any information on the external acceleration, i.e., . 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 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
As and are of the same order, the pseudoelements should then satisfy
Then, the total velocity could be expressed as and the osculating velocity and the position could be expressed as
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:
2.4.3. Pseudoequinoctial Fourier Series Method
Eventually, within the transfer time , the constraint and velocity increment can be separately expressed as
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 where denote the relaxation factors. The relationship between the fuel cost and final velocity increment is expressed as [18] where denotes the specific impulse and 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 PEFFS 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. TwoBody Rendezvous Trajectory Design
In this section, the twobody rendezvous trajectory design using PEFFS is investigated. It is assumed that in transfer problems including low Earth orbit (LEO), EarthMars, Earthasteroid, and Earthcomet, 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 nearEarth space, the distance unit (DU) refers to the average radius of the Earth, and the time unit (TU) refers to the period of a circular orbit with the radius of the Earth. For the deepspace exploration problems, DU refers to the astronomical unit (AU), i.e., . TU equals the period of an orbit with the radius of 1 DU around the Sun, i.e., . 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 i76700 processor, a Windows 7 operating system, and 8 GB RAM was used for the simulation.
The data of the deepspace 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 Tempel1 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 , , were set as 10, 1, and 1, respectively, by trial and error; and the error tolerance was .

During the calculation, the revolution number in 3DFFS is traversed. Some initial determination was introduced to narrow the search range of the number of revolution. The relationship between the maximum acceleration , phase angle , and transfer time , was considered. The phase angle can be expressed as 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 needs to be added. 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
At the same time, considering the thruster capability, the minimum number of revolutions is where and denote the velocities of the initial and final trajectories, and and denote the semimajor axes of the initial and final trajectories, respectively. The physically impossible condition 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, and can be rounded instead of the upper and lower roundedoff 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 of ; 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.

The 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 longterm 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.
The above comparison indicates that the proposed method has the lowest acceleration peak of . 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 PEFFS and 3DFFS 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 of 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.

Because the number of terms is adjusted in the shape library, the changes in PEFFS are smoother, and there are more feasible solutions in the whole process than those in 3DFFS. The optimal results of the trajectory and thrust acceleration of the transfer are shown in Figures 8 and 9, respectively.
The best launch window is 7955 MJD2000 with 440 days TOF and one complete revolution around the Sun. The corresponding shape number in Table 4 is five, and the peak value of the thrust acceleration is . The comparison with other methods is presented in Table 5.
The average computational time of PEFFS and 3DFFS 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 3DFFS 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 PEFFS has a broader convergence field and provides more feasible solutions.
3.3. Earth to Comet Tempel1
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. is set as , 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.

Figure 10 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.
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 is 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.
The velocity increment of PEFFS is better than that of 3DFFS. At the same time, it is ensured that the thrust acceleration satisfies the constraint conditions. The average computational time of PEFFS and 3DFFS 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 shapebased method is suitable for a largedesign space research.
From the analysis of the above nearEarth trajectory design and deepspace exploration trajectory design, it can be concluded that PEFFS can be applied to design rendezvous trajectories from nearEarth to deep space, and the obtained thrust acceleration changes are smoother and slower than those obtained by 3DFFS. 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 3DFFS, with a gradual increase in the eccentricity and inclination angle of the transfer trajectory, more feasible solutions can be obtained than those obtained by 3DFFS. For trajectory rendezvous problems with large eccentricity and inclination angle, improved velocity increment results can be obtained using PEFFS. In addition, the introduction of a different number of terms in Tables 2, 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 3DFFS. 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. LowThrust Trajectory Design and Analysis for CRTBP
In the above cases, transfer processes are frequently affected by the central gravitational body; hence, the twobody 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 EarthMoon rendezvous trajectory design as an example, the CRTBP, a method that employs two twobody models to calculate the threebody lowthrust trajectory, is proposed. The free intermediate section without lowthrust acceleration during the EarthMoon 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 EarthMoon transfer was designed using PEFFS and compared to the 3DFFS 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] where denotes the EarthMoon distance, and and 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 is the angle between the velocity vector and the position vector at capture point , where . The position and velocity vectors at perilune are denoted as and , 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 degenerates to .
When a Kepler orbit is considered, the semilatus rectum of the trajectory can be calculated using the critical capture parameters and . The semimajor axis , eccentricity , and phase angle are expressed as [37] where denotes the energy factor, which can be expressed as where denotes the lunar gravitational constant. The perilune radius can be obtained from equation (43). When 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 to point is where
Except for the case of a hard landing, two capture criteria are given:
Criterion 1. During the transfer process, the range of feasible radius 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
When the initial state of the vehicle is propagated without external thrust, the relationship between the capture parameters and and perilune distance , and that between and and transfer time 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.
The above figures indicate that if is within the range of 0.1 to 0.1 or 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 is larger, the vehicle will fly out of the LSOI during the transfer time . As increases, the feasible range of shrinks gradually. Within the range of larger than 385 m/s, the and transfer time can be feasible only under the condition of an extremely small , but the vehicle does not satisfy Criterion 2 and will fly out of the LSOI after .
To expand the area satisfying Criterion 2, thrust is applied to brake the vehicle. Constant thrust acceleration from to is imposed opposite the direction of the velocity within a time of . 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.
From Figures 16 and 17, it can be inferred that the feasible range enlarges due to the deceleration thrust. The darkcolored area indicates that the vehicle is automatically captured. The lightcolored 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 , the thrust for the lunar capture has to increase. decreases with smaller . In the area between 2 and 1, is at the edge of the LSOI and the amplitude of deceleration is large.
4.2. EarthMoon Transfer Launch Window Design
The geocentric twobody model was adopted in the escape section, and the proposed PEFFS 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
During the transfer process, the trajectory will pass through the LSOI. At that moment, the position and velocity relative to the Moon can be expressed as
The angle between the position and the velocity can be expressed as
The range of the angle is between and .
The first step is to design the whole trajectory with the initial GEO states and the final motion states of the Moon by PEFFS under the twobody dynamics. Because the treebody 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 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 . 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. refers to the velocity increment from the GEO to the LSOI. refers to the total velocity increment under the twobody dynamics. The proportion is the ratio of the to the .

Within 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.
The 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 2DFFS, due to the deficiency of the intermediate phase. However, the computational time of PEFFS 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 lowthrust 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 PEFFS are shown in Figures 18 and 19, respectively.
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/s^{2}. 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 EarthMoon 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 PEFFS can still be applied.
5. Conclusions
In this study, a new method for the shapebased trajectory design was proposed. The inverse dynamic model for the twobody system was established and the relationship between acceleration and the pseudoequinoctial elements was proposed. Following this, the proposed method was examined using twobody system problems for deepspace explorations such as EarthMars, Earth1989ML, and EarthTempel1. 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 threebody trajectory transfer problems, the proposed method combined with the patched conics method was presented to analyze the lowthrust EarthMoon transfer window. The combined method was applied to comprehensively design the initial lowthrust geocentric trajectory for EarthMoon transfer, which can improve the design efficiency and ensure precision at the same time.
Nomenclature
Greek Symbols:  Energy factor 
:  Capture angle (rad) 
:  Phase angle (rad) 
:  Elements of the orbit motion states 
:  Gravitational constant of the central body (m^{3}/s^{2}) 
:  Argument of perigee (rad) 
:  Right ascension of ascending node (rad). 
:  Semimajor axis (m) 
c_{1}:  Relaxation factor of the acceleration constraints 
c_{2}:  Relaxation factor of the angular velocity constraints 
c_{3}:  Relaxation factor of the trajectory direction constraint 
:  Eccentricity 
:  The second modified equinoctial element (rad) 
:  acceleration vector (m/s^{2}) 
:  state transform function 
:  Thrust acceleration vector (m/s^{2}) 
:  The third modified equinoctial element (rad) 
:  Gravitational acceleration (m/s^{2}) 
:  Gravity loss 
:  Control transform function 
:  The fourth modified equinoctial element (rad^{2}) 
:  Observation transform function 
:  Inclination (rad) 
:  Specific impulse (s) 
:  Jacobian matrix 
:  The fifth modified equinoctial element (rad^{2}) 
:  Coefficients of each sine function term 
:  Coefficients of each cosine function term 
:  The last modified equinoctial element (rad) 
:  Lie derivative function of 
:  Lie derivative function of 
:  Mass (kg) 
:  Number of terms 
:  The first modified equinoctial element (m) 
:  Radius vector (m) 
:  Time of flight (s) 
:  Control vector (m/s^{2}) 
:  Velocity vector (m/s) 
:  Velocity increment (m/s) 
:  State variables 
:  Measured variables 
:  axial distance (m). 
3D:  Three dimensional 
AU:  Astronomical unit 
CRTBP:  Circular restricted threebody problem 
DU:  Distance unit 
ESA:  European Space Agency 
FFS:  Finite Fourier series 
GEO:  Geosynchronous orbit 
GPOPS:  General pseudospectral optimization software 
GTOC 2:  The 2^{nd} Global Trajectory Optimization Competition 
IP:  Inverse polynomial method 
JAXA:  Japan Aerospace Exploration Agency 
JPL:  Jet Propulsion Laboratory 
LEO:  Low Earth orbit 
LSOI:  Lunar sphere of influence 
MEEs:  Modified equinoctial elements 
MJD:  Modified Julian date 
NASA:  National Aeronautics and Space Administration 
PE:  Pseudoequinoctial 
SQP:  Sequential quadratic programming 
TOF:  Time of flight 
TPBV:  Twopoint boundary value 
TU:  Time unit. 
:  Capture point 
:  Perilune 
:  Desired variable 
:  Earth 
:  Earth to Moon 
error:  Deviation of the variable 
escape:  Escape trajectory 
:  Final variable 
gau:  Gauge variable 
impulse:  Impulse trajectory 
:  Variable relative to Moon 
lowthrust:  Lowthrust trajectory 
:  Moon 
max:  Maximum value 
min:  Minimum value 
:  Number 
:  Normal component in RTN coordinates 
:  Initial variable 
osc:  Osculating variable 
:  Radial component in cylindrical coordinates 
rev:  Revolution 
tof:  Time of flight 
:  Axis component in inertial coordinates 
:  Axis component in inertial coordinates 
:  Vertical component in cylindrical coordinates/axis component in inertial coordinates 
:  Tangential component in cylindrical coordinates. 
~:  Pseudoequinoctial elements 
:  Firstorder derivative with respect to time 
:  Secondorder derivative with respect to time 
:  Variable relative to Moon. 
Data Availability
The data used to support the findings of this study (No. 1364834) may be released upon application to the College of Aerospace Science and Engineering, National University of Defense Technology, which can be contacted at lihaiyang@nudt.edu.cn.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
We would like to thank the anonymous reviewers for their invaluable feedback. This study was cosupported by the National Natural Science Foundation of China (Nos. 11472301 and 11372345).
References
 M. D. Rayman and S. N. Williams, “Design of the first interplanetary solar electric propulsion mission,” Journal of Spacecraft and Rockets, vol. 39, no. 4, pp. 589–595, 2002. View at: Publisher Site  Google Scholar
 C. T. Russel, F. Capaccioni, A. Corradini et al., “Dawn discovery mission to Vesta and Ceres: present status,” Advances in Space Research, vol. 38, no. 9, pp. 2043–2048, 2006. View at: Publisher Site  Google Scholar
 T. B. John and O. E. Sven, “Optimal low thrust trajectories to the moon,” SIAM Journal of Applied Dynamical Systems, vol. 2, no. 2, pp. 144–170, 2003. View at: Publisher Site  Google Scholar
 Y. Langevin, “Chemical and solar electric propulsion options for a cornerstone mission to mercury,” Acta Astronautica, vol. 47, no. 2, pp. 443–452, 2000. View at: Publisher Site  Google Scholar
 “2^{nd} Global Trajectory Optimisation Competition (GTOC2),” Nov, 2008, https://sophia.estec.esa.int/gtoc_portal/?page_id=15. View at: Google Scholar
 P. Enright and B. Conway, “Optimal finitethrust spacecraft trajectories using collocation and nonlinear programming,” Journal of Guidance, Control, and Dynamics, vol. 14, no. 5, pp. 981–985, 1991. View at: Publisher Site  Google Scholar
 J. L. Junkins and T. Ehsan, “Exploration of alternative state vector choices for lowthrust trajectory optimization,” Journal of Guidance, Control, and Dynamics, vol. 42, no. 1, pp. 47–64, 2018. View at: Publisher Site  Google Scholar
 J. T. Betts, “Survey of numerical methods for trajectory optimization,” Journal of Guidance, Control, and Dynamics, vol. 21, no. 2, pp. 193–207, 1998. View at: Publisher Site  Google Scholar
 P. R. Ryan, “Primer vector theory applied to global lowthrust trade studies,” Journal of Guidance, Control, and Dynamics, vol. 30, no. 2, pp. 460–472, 2007. View at: Publisher Site  Google Scholar
 A. L. Herman and B. A. Conway, “Direct optimization using collocation based on highorder GaussLobatto quadrature rules,” Journal of Guidance, Control, and Dynamics, vol. 19, no. 3, pp. 592–599, 1996. View at: Publisher Site  Google Scholar
 A. V. Rao, D. A. Benson, M. A. P. Christopher Darby, C. Francolin, I. Sanders, and G. T. Huntington, “Algorithm 902: GPOPS, a MATLAB software for solving multiplephase optimal control problems using the gauss pseudo spectral method,” ACM Transactions on Mathematical Software, vol. 37, no. 2, article 22, 2010. View at: Publisher Site  Google Scholar
 H. S. Tsien, “Takeoff from satellite trajectory,” Journal of the American Society, vol. 23, no. 4, pp. 233–236, 1953. View at: Google Scholar
 D. J. Benny, “Escape from a circular orbit using tangential thrust,” Journal of Jet Propulsion, vol. 28, no. 3, pp. 167–169, 1958. View at: Publisher Site  Google Scholar
 D. F. Lawden, “Optimal escape from a circular trajectory,” Acta Astronautica, vol. 4, pp. 218–233, 1958. View at: Google Scholar
 A. E. Petropoulos and J. M. Longuski, “Shapebased algorithm for the automated design of lowthrust, gravity assist trajectories,” Journal of Spacecraft and Rockets, vol. 41, no. 5, pp. 87–796, 2004. View at: Publisher Site  Google Scholar
 B. J. Wall and B. A. Conway, “Shapebased approach to lowthrust rendezvous trajectory design,” Journal of Guidance, Control, and Dynamics, vol. 32, no. 1, pp. 95–101, 2009. View at: Publisher Site  Google Scholar
 A. E. Petropoulos and J. A. Sims, “A review of some exact solutions to the planar equations of motions of a thrusting spacecraft,” in Proceedings of the 2nd International Symposium on LowThrust Trajectory, Toulouse, France, 2002. View at: Google Scholar
 B. J. Wall, “Shapebased approximation method for lowthrust trajectory optimization,” in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, Hawaii, 2008. View at: Publisher Site  Google Scholar
 P. De and M. Vasile, “Preliminary design of lowthrust multiple gravityassist trajectories,” Journal of Spacecraft and Rockets, vol. 43, no. 5, pp. 1065–1076, 2006. View at: Publisher Site  Google Scholar
 D. M. Novak and M. Vasile, “Improved shaping approach to the preliminary design of lowthrust trajectories,” Journal of Guidance, Control, and Dynamics, vol. 34, no. 1, pp. 128–147, 2011. View at: Publisher Site  Google Scholar
 E. Taheri, K. Ilya, and A. Ella, “Shaping lowthrust trajectories with thrusthandling feature,” Advances in Space Research, vol. 61, no. 3, pp. 879–890, 2018. View at: Publisher Site  Google Scholar
 D. J. Gondelach and R. Noomen, “Hodographicshaping method for lowthrust interplanetary trajectory design,” Journal of Spacecraft and Rockets, vol. 52, no. 3, pp. 728–738, 2015. View at: Publisher Site  Google Scholar
 M. Razzaghi, “Optimal control of linear timevarying systems via Fourier series,” Journal of Optimization Theory and Applications, vol. 65, pp. 375–384, 1990. View at: Publisher Site  Google Scholar
 R. V. Dooren and J. Vlassenbroeck, “A Chebyshev technique for solving nonlinear optimal control problems,” IEEE Transactions on Automatic Control, vol. 33, pp. 333–339, 1988. View at: Publisher Site  Google Scholar
 M. Razzaghi, “Fourier series direct method for variational problems,” International Journal of Control, vol. 48, pp. 887–895, 1988. View at: Publisher Site  Google Scholar
 Q. Fang, X. Wang, C. Sun, and J. Yuan, “A shapebased method for continuous lowthrust trajectory design between circular coplanar orbits,” International Journal of Aerospace Engineering, vol. 2017, Article ID 9234905, 10 pages, 2017. View at: Publisher Site  Google Scholar
 M. Y. Huo, G. Zhang, N. Qi, Y. Liu, and X. Shi, “Initial trajectory design of electric solar wind sail based on finite Fourier series shapebased method,” IEEE Transactions on Aerospace and Electronic Systems, p. 1, 2019. View at: Publisher Site  Google Scholar
 O. Abdelkhalik and E. Taheri, “Approximate onoff lowthrust space trajectories using Fourier series,” Journal of Spacecraft and Rockets, vol. 49, no. 5, pp. 962–965, 2012. View at: Publisher Site  Google Scholar
 E. Taheri and O. Abdelkhalik, “Fast initial trajectory design for lowthrust restrictedthreebody problems,” Journal of Guidance, Control, and Dynamics, vol. 38, no. 11, pp. 2146–2160, 2015. View at: Publisher Site  Google Scholar
 E. Taheri and O. Abdelkhalik, “Shape based approximation of constrained lowthrust space trajectories using Fourier series,” Journal of Spacecraft and Rockets, vol. 49, no. 3, pp. 535–545, 2012. View at: Publisher Site  Google Scholar
 E. Taheri and O. Abdelkhalik, “Initial threedimensional lowthrust trajectory design,” Advances in Space Research, vol. 57, pp. 890–903, 2016. View at: Publisher Site  Google Scholar
 J. T. Betts, “Very lowthrust trajectory optimization using a direct SQP method,” Journal of Computational and Applied Mathematics, vol. 120, no. 1, pp. 27–40, 2000. View at: Publisher Site  Google Scholar
 Y. Gao, Advances in LowThrust Trajectory Optimization and Flight Mechanics, University of Missouri Columbia, Columbia, 2003.
 A. Giulio, T. Douglas, and T. Alberto, “Model predictive control architecture for rotorcraft inverse simulation,” Journal of Guidance, Control, and Dynamics, vol. 36, no. 1, pp. 207–217, 2013. View at: Publisher Site  Google Scholar
 L. Linghai, Inverse Modelling and Inverse Simulation for System Engineering and Control Applications, Department of Electrical Engineering, University of Glasgow, 2007.
 “Jet Propulsion Laboratory (JPL) ephemeris,” https://ssd.jpl.nasa.gov. View at: Google Scholar
 A. V. David, Fundamentals of Astrodynamics and Applications, Microcosm Press, EI Segundo, CA, USA, 2nd edition, 2004.
 B. L. Pierson and C. A. Kluever, “Threestage approach to optimal lowthrust earthmoon trajectories,” Journal of Guidance, Control, and Dynamics, vol. 17, no. 6, pp. 1275–1282, 1994. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Wanmeng Zhou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.