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.
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 -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 -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 where refers to the gravitational constant of the central body and refers to the acceleration produced by the external force. In traditional shape-based 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 two-body dynamic equation is expressed by the modified equinoctial elements in the Gaussian perturbation equation [19]: where , , and 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 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 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:
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 first-order and second-order 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. 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 where refers to the observation function. According to the current states obtained in the real-time 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 third-order polynomial that satisfies the boundary conditions:
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. 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 Moore-Penrose 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 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 the 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., . 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 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 , , were set as 10, 1, and 1, respectively, by trial and error; and the error tolerance was .
During 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 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 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 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 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.
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 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 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 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 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 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. 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 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 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 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] where denotes the Earth-Moon 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 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 , 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. 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
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 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 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 two-body 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 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.
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.
Nomenclature
Greek Symbols| : | Energy factor |
| : | Capture angle (rad) |
| : | Phase angle (rad) |
| : | Elements of the orbit motion states |
| : | Gravitational constant of the central body (m3/s2) |
| : | Argument of perigee (rad) |
| : | Right ascension of ascending node (rad). |
| : | Semimajor axis (m) |
| c1: | Relaxation factor of the acceleration constraints |
| c2: | Relaxation factor of the angular velocity constraints |
| c3: | Relaxation factor of the trajectory direction constraint |
| : | Eccentricity |
| : | The second modified equinoctial element (rad) |
| : | acceleration vector (m/s2) |
| : | state transform function |
| : | Thrust acceleration vector (m/s2) |
| : | The third modified equinoctial element (rad) |
| : | Gravitational acceleration (m/s2) |
| : | Gravity loss |
| : | Control transform function |
| : | The fourth modified equinoctial element (rad2) |
| : | Observation transform function |
| : | Inclination (rad) |
| : | Specific impulse (s) |
| : | Jacobian matrix |
| : | The fifth modified equinoctial element (rad2) |
| : | 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/s2) |
| : | Velocity vector (m/s) |
| : | Velocity increment (m/s) |
| : | State variables |
| : | Measured variables |
| : | axial distance (m). |
| 3D: | Three dimensional |
| AU: | Astronomical unit |
| CRTBP: | Circular restricted three-body problem |
| DU: | Distance unit |
| ESA: | European Space Agency |
| FFS: | Finite Fourier series |
| GEO: | Geosynchronous orbit |
| GPOPS: | General pseudospectral optimization software |
| GTOC 2: | The 2nd 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: | Two-point 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: | Low-thrust 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 |
| : | First-order derivative with respect to time |
| : | Second-order 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 [email protected].
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).