Abstract
The shape-based method can provide suitable initial guesses for trajectory optimization, which are useful for quickly converging a more accurate trajectory. Combined with the optimal control theory, an optimized shape-based method using the finite Fourier series is proposed in this paper. Taking the flight time-fixed case and the time-free case into account, respectively, the optimized shape-based method, which considers the first-order optimal necessary conditions, can guarantee that not only an orbit designed during the preliminary phase is optimal, but also the thrust direction is not constrained to be tangential. Besides, the traditional shape-based method using the finite Fourier series, in which the thrust direction is constrained to be tangential, is developed for the time-free case in this paper. The Earth-Mars case and the LEO-GEO case are used to verify the optimized shape-based method’s feasibility for time-fixed and time-free continuous low-thrust trajectory design between circular coplanar orbits, respectively. The optimized shaped-based method can design a lower cost trajectory.
1. Introduction
Recently, continuous low-thrust trajectory design and optimization are becoming increasingly popular [1, 2], although they are very challenging and time-consuming. What is particular is that the continuous low-thrust trajectory design consists of two phases: preliminary design and precise design [3]. To pursue a faster optimization and a more accurate trajectory, the preliminary design phase is expected to provide an efficient initial guess for trajectory optimizers. The shape-based (SB) method is one of the most efficient methods during this preliminary design. The SB method assumes that some functions contain a spacecraft’s trajectory, and therefore boundary conditions are used to calculate the parameters of the functions, thus analytically obtaining the needed thrust during the spacecraft’s flight.
Many kinds of SB methods have been proposed by researchers. For instance, Petropoulos and Longuski [4, 5] developed an exponential sinusoid (ES) method for the two-dimensional (2D) interplanetary transfer trajectory design. The ES method constrains the thrust direction to be tangential. Izzo [6] utilized this method to investigate the multirevolution Lambert’s problem and simplified the interplanetary low-thrust trajectory design procedure. Cui et al. [7] proposed a new search approach algorithm for the launch window of low-thrust gravity-assist missions based on the ES method, which has fewer searching variables and is more efficient than the traditional SB methods. But the ES method cannot satisfy the full consideration conditions of circular terminal orbits unless thrusters provide impulsive propulsion; besides, the parameters of a shape cannot be solved when other constraints are introduced.
Zheng et al. [8] proposed a new trajectory shape called the logarithmic spiral-based (LS) non-Keplerian orbit. The feasibility and essential characteristics of the LS non-Keplerian orbit are analyzed. The analytical geocentric distance expression and the phase angle expression about flight time subject to a tangential thrust are derived. But the LS method cannot satisfy the terminal constraints as well as the ES method.
To overcome the above-mentioned disadvantages, Wall and Conway [9, 10] developed a 2D inverse polynomial (IP) method. The fifth-order IP method can be used to design the transfer trajectory for the time-free case, while the sixth-order IP method is designed for the time-fixed case. But the thrust direction is also constrained to be tangential in the IP method, which cannot handle the thrust constraint very well. Shang et al. [11] proposed a semianalytical Lambert algorithm based on the -degree IP method in order to improve the precision of preliminary design for an interplanetary low-thrust transfer trajectory. Considering thrust and radius constraints, Wang et al. [12] proposed a modified IP method for both the time-free transfer case and the time-fixed rendezvous case. Compared with the original IP method, the modified one can satisfy the thrust and radius constraints through optimizing the polynomial orders. To realize low-thrust trajectory design between elliptical orbits, Xie et al. [13] constructed a new shape function about two semimajor-axis parameters, which are polynomials in the polar angle. With tangential thrust, a fifth-order and sixth-order method is designed for time-free and time-fixed cases as in the IP method.
Taheri and Abdelkhalik [14] and Abdelkhalik and Taheri [15] proposed a new shape-based trajectory design method using the finite Fourier series. With the hypothesis of tangential thrust, a preliminary trajectory that satisfies the maximum thrust constraint is designed with this SB trajectory design method.
Shaping the velocity components, Gondelach et al. [16] proposed a novel low-thrust trajectory design method called hodographic-shaping (HS) method. These velocity functions are assumed to be some sets of simple base functions. Extra parameters are used to make the trajectory design and optimization more flexible.
In a word, all recent SB methods (except the HS method) design spacecraft trajectories based on the tangential thrust assumption and cannot guarantee that the designed trajectory is optimal without the first-order optimal necessary conditions. Meanwhile, almost SB methods require iterative calculations or constraint optimization to match the total flight time. However, it is difficult to determine a suitable flight time during or before the preliminary design phase.
Therefore, combined with optimal control theory, an optimized shape-based method using finite Fourier series, which can easily overcome the above-mentioned shortcomings, is proposed in this paper. Regarding spacecraft three-dimensional (3D) trajectory design, Wall [9], Novak and Vasile [17], and Taheri [18] presented their study advances. However, 3D trajectory design is not the key point in this paper, because the 2D case is enough to illustrate the idea of our method.
The paper is organized as follows. In Section 2, the spacecraft dynamics model in polar coordinate is developed. Then, the proposed method is introduced in Section 3. In the time-fixed rendezvous case and time-free transfer case, respectively, the first-order necessary conditions are derived from the Hamiltonian function. Through expanding the state variables into expressions of finite Fourier series, the optimal control problem is converted to a nonlinear programming (NLP) problem that involves Fourier series coefficients. The process of the proposed method is given in Section 4. In Section 5, two examples are used to verify the optimized shape-based method’s feasibility for continuous low-thrust trajectory design between circular coplanar orbits, and the advantage that the proposed method can design a lower cost trajectory is proven by comparing it with other SB methods.
2. Spacecraft Orbital Model
In polar coordinate, the spacecraft’s orbital motion model without considering any perturbation and celestial bodies’ rotation is established as where superscript “” indicates a derivative with respect to time ; is the magnitude of the position vector; is the polar angle; is the magnitude of spacecraft radial velocity vector; is the magnitude of its circumferential velocity vector; is the gravitational parameter; is the thrust acceleration magnitude of the spacecraft; is its steering angle, and is its flight path angle, as shown in Figure 1.
Instead of shaping variables as a function of time, the polar angle also can be used as an independent variable. In that case, the variables have to be analytically integrable over to obtain the change in position.
In this paper, a new thrust acceleration parameter is introduced, its direction is the same as , and its magnitude is defined as shown in
Instead of derivatives with respect to time , these state variables themselves are derivatives with respect to the polar angle . According to the second equation of (1), the spacecraft’s orbital motion model can be reestablished aswhere suffix “” indicates a derivative with respect to polar angle .
3. The Optimized Shape-Based Method Using Fourier Series
3.1. Performance Index and Boundary Conditions
Usually the performance index for the low-thrust trajectory design is set to minimize the flight time or fuel consumption. In this paper, we consider the minimal characteristic velocity (i.e., minimal fuel consumption), as shown in
To accomplish a spacecraft’s transfer successfully, some constraints such as those given in (5) should be satisfied.where represent the spacecraft’s position, velocity, and acceleration conditions, respectively; represent initial and terminal condition, respectively.
Actually, the fact that mentioned in (5) are regarded as generalized boundary conditions, which may be the first-order or second-order derivative of a certain variable, is more reasonable.
3.2. The Flight Time-Fixed Case
3.2.1. Traditional Fourier Series (TFS) Method
In the time-fixed case, it is assumed that the thrust is aligned along or against the velocity vector; that is, , where . The TFS method for the time-fixed case was studied in detail in Taheri’s doctoral dissertation [18]. Besides, we only consider the unconstrained version of the finite Fourier series method in this paper; because the thrust magnitude constraint is not a key point, we focus on optimizing the FS method here.
The following equation is derived from the fourth equation of (1):
Substituting (6) into the third equation of (1), the following equation is derived:where the tangential thrust assumption can be written as
Substituting the first and second equations of (1) and the tangential thrust assumption into (7), one can be rewritten as
According to Fourier series expansion, the radius and the polar angle can be approximated as follows:where means ; is the number of finite Fourier terms; are Fourier coefficients; is the total flight time.
Substituting the state approximation (10) into (9), the differential equation is converted to a nonlinear algebraic equation, in which the only unknowns are the Fourier coefficients and the independent time variable:
This shows that, in the flight time-fixed case that considers the tangential thrust direction hypothesis, (4), (5), and (11) become an NLP problem about Fourier coefficients.
3.2.2. Optimized Fourier Series (OFS) Method
In light of (1) and (4), the Hamilton function is written as
Based on the optimal control theory, costate equations and control equations are shown as
Let and ; the third and fourth equations of (1) are rewritten as
With (15), the spacecraft’s thrust acceleration magnitude is solved.
Meanwhile, the thrust acceleration direction is defined by the results of and .
The following equation is derived from (14) and (15):
The expressions of the costate variables are obtained with the solution of (17):
According to the second and fourth equations of (13), the value of constant can be represented as
The following equation is rearranged from the first equation of (13):where the costate variable is solved with the third equation of (13)
Therefore, all the costate variables can be expressed as the functions of state variables .
Substituting the state approximations of (10) into (20), the differential equation is converted to a nonlinear algebraic equation, in which the only unknowns are the Fourier coefficients and the independent time variable:
This shows that, in the flight time-fixed with no limitation on thrust direction, (4), (5) and (22) become an NLP problem about Fourier coefficients.
3.3. The Flight Time-Free Case
3.3.1. TFS Method
In the flight time-free case, it is assumed that the thrust is aligned along or against the velocity, the same as what was mentioned in Section 3.2.1. This method for the flight time-free case is a new study in this paper.
The following equation is derived from the second and third equations of (3):where the tangential thrust assumption is expressed as
Substituting the first equation of (3) and (24) into (23), the following equation is derived:
Dividing (25) by , the following equation is derived:
Considering the thrust direction assumption, (26) is rewritten as
Equation (27) can be simplified as
According to the Fourier series expansion, the radius and the circumferential velocity magnitude can be approximated as follows:where means ; ; is the number of revolutions around the attracting central body, as shown in Figure 2.
Substituting the state approximations of (29) into (28), the differential equation is converted to a nonlinear algebraic equation, in which the only unknowns are the Fourier coefficients and the polar angle:
This shows that, in the flight time-free case that considers the tangential thrust direction hypothesis, (4), (5), and (11) become an NLP problem about Fourier coefficients.
According to (3), is represented aswhere
The spacecraft’s flight time is solved with the following equation:
3.3.2. OFS Method
In light of (3) and (4), the Hamiltonian function is expressed as
Based on the optimal control theory, the costate equations and control equations are expressed as
With , (35) and (3) can be rewritten as
Let and ; (38) is rewritten as
The magnitude of is solved with (39):
Meanwhile, the direction of is defined by the results of and .
The following equations are derived from (36) and (39):
The expressions of the costate variables are obtained with the solution of (41):
The following equation is derived from the first and third equations of (37):where the costate variable is solved with the second equation of (37) and (42):
So all the costate variables can be expressed as functions of state variables .
Substituting the state approximation of (29) into (43), the differential equation is converted to a nonlinear algebraic equation, in which the only unknowns are the Fourier coefficients and the polar angle:
This shows that, in the flight time-free case with no limitation on thrust direction, (4), (5), and (45) become an NLP problem about Fourier coefficients.
4. Transfer Orbit Design with the FS Method
The flowchart of the spacecraft’s transfer trajectory design with the FS method is shown in Figure 3, in which FCs is short for “Fourier coefficients.”
The fmincon function is a MATLAB command used to solve multivariable, nonlinear and optimal problem with constraints. However, it requires that a set of initial guesses for the Fourier coefficients should be obtained. Reference [14] proposed some rough approaches used to gain the initial guesses. Using the FS method for designing an interplanetary low-thrust trajectory, a cubic polynomial function is used to obtain initial guesses for the Fourier coefficients. The constraint about finite Fourier series terms is set to satisfy the boundary conditions. Although there is no upper limit on the number of included Fourier terms, the computational efficiency and precision are important considerations.
Based on the initial guesses, the new Fourier coefficients are calculated with the fmincon function. Then, the analytical trajectory and the thrust expressed with the finite Fourier coefficients are found. These expressions are used to offer initial guesses for detailed trajectory optimizers. After the preliminary design phase, trajectory optimization is required to show the advantage of the efficient initial guess from OFS method compared with the guess from TFS.
In this paper, we use the direct collocation method [19]. In the direct collocation method, the optimization model, performance index, and constraints are the same as those in other SB methods, for example, (1) or (3), (4), and (5). The optimal control problem would be converted into an NLP problem. The total flight interval is discretized into intervals, and two endpoints of each interval are called “node.” In this paper, the selection of the number of nodes is not a key point. Certainly, the simulation results will be more accurate as the number of nodes increases.
5. Simulation Examples
In the time-fixed and time-free cases, the continuous low-thrust Earth-Mars rendezvous and LEO-GEO transfer are studied, respectively. The simulation of these cases has been performed with MATLAB 2014a on an Intel Core i5 2.6 GHz computer with Windows 8.
5.1. The Earth-Mars Rendezvous
The FS method for the time-fixed case can be applied to the interplanetary exploration, asteroid deflection, rendezvous and docking missions, and so on. In the time-fixed case, we design the continuous low-thrust Earth-Mars rendezvous trajectory by using canonical units, where 1 distance unit (DU) is 1 AU and time unit (TU) is 1 year. The boundary conditions and input parameters are listed in Table 1, in which is the number of nodes in the direct collocation method.
Figure 4 gives the spacecraft’s rendezvous orbits using different FS methods. Figure 5 shows the angle relation between velocity direction and acceleration direction using the OFS method.
Different from other SB methods, which typically suppose that the acceleration direction of a spacecraft is parallel with its velocity direction, the OFS method does not need to constrain the relation between acceleration direction and velocity direction as shown in Figures 4 and 5.
In Figure 5, means the difference between steering angle and flight path angle; for example, . The figure shows that the difference fluctuates within the range of −2° to 2° during the large part of flight time, and the maximum difference is −12.59°. The acceleration direction varies sharply near the initial and terminal points, because the initial and terminal conditions result in the obvious change of the two parameters and .
Figure 6 gives the spacecraft’s thrust acceleration profiles using different FS methods. In this case, we can observe that the thrust acceleration profile using the OPS method is smoother than the thrust acceleration profile using the TFS method.
Regarding the thrust direction, the proposed method does not assume it to be tangential. When the obtained solution is used as initial guess in an optimizer, Figure 7 shows the thrust direction after optimization in the Earth-Mars case. The figure shows that the difference fluctuates within the range of −4° to 4° during the large part of flight time, and the maximum difference is −15°. Similar to the preliminary phase, the thrust direction after optimization varies sharply near the initial and terminal points. It justify that it is significant to relax the thrust direction or constrain it.
Table 2 gives simulation results using different initial guesses. It shows evidently that the OFS method obtains a lower cost rendezvous trajectory, and this verifies the applicability of our OFS method to continuous low-thrust trajectory design. Whether during the preliminary design phase or the precise design phase, the transfer trajectory designed with the OFS method consumes minimum fuel, but the optimization time of the OFS method is a little longer.
5.2. The LEO-GEO Transfer
The FS method for the time-free case can also be applied for a spacecraft orbit raise and maneuver mission. For the time-free case, the LEO-GEO transfer mission is considered, where 1 DU is 1 and 1 TU is 808.67 s. In this case, the flight time is unknown to us previously, but the terminal position is given. Table 3 lists the boundary conditions and input parameters.
In the time-free case, and their first-order derivatives with respect to polar angle shall be calculated as new boundary conditions, as shown in
Figure 8 gives the spacecraft’s LEO-GEO transfer orbits using different FS methods. Figure 9 shows the angle relation between velocity direction and acceleration direction using the OFS method. Figure 10 gives the spacecraft’s thrust acceleration profiles using different FS methods. Figure 11 shows the thrust direction after optimization in the LEO-GEO case.
Figures 8 and 10 prove that the shape-based trajectory design method using the finite Fourier series can be used in the flight time-free case. This paper extends the scope of application of the Fourier series method to the flight time-free case. Figure 10 shows that, in the flight time-free case, different FS methods have the similar trend: during the initial phase of the flight, the thrust acceleration varies obviously and then tends to be stable.
Figures 8 and 9 also show that, similar to the Earth-Mars transfer in the flight time-fixed case, the OFS method does not need to constrain the angle relation between acceleration direction and velocity direction. As shown in Figure 9, the phenomena that the acceleration direction varies sharply near the initial and terminal points also exist, because the parameters are and at the initial point and terminal point, respectively, which result in and . As shown in Figures 5 and 9, it is obvious to find that the phenomena where the thrust direction is not tangential near the initial and terminal points are determined by the OFS method itself, because of and . As shown in Figure 11, some similar results, which are shown in Earth-Mars rendezvous case, can be found.
Table 4 presents simulation results using different initial guesses. As shown in Table 4, with the OFS method, the spacecraft has the longest flight time, and with the TFS method, its flight time is almost the same as that with the OFS method because the rough approach mentioned in Section 4, which can offer initial Fourier coefficients, leads to the approximate trajectory of the spacecraft. With the rough approach, the flight time is 112.2044 TU, which is almost the same as that with the OFS method and TFS method. Therefore, it suggests that the rough approach, which offers initial coefficients for the FS method, plays an important role.
The same as in the time-fixed case, the OFS method in the time-free case can provide better initial guesses for optimizers, and the transfer trajectory designed with the OFS method consumes minimal fuel, which is 0.5880 DU/TU before optimization and 0.5803 DU/TU after optimization. For all FS methods, the flight time after optimization becomes slightly longer. The increment with the OFS method is 0.0811 TU and TFS method is 0.2371 TU.
With regard to the optimization time, we obtain the similar simulation results that the flight time under the OFS method is a little longer than that under the TFS method, as shown in Tables 2 and 4. From the derivation in Section 3, we know that a more complicated NLP problem of the OFS method needs to be solved than the NLP problem of TFS method.
5.3. Remark
In this subsection, two important points are given to illustrate the OFS method.
With regard to the number of terms in the Fourier expansion, theoretically, the finite Fourier series expression will be very close to the real function as the number of Fourier series terms increases according to the characteristic of series theory. So the OFS method can get analytical trajectory and thrust which is very close to the optimal solution, if the number of series terms is enough. Actually, in this case, OFS looks more like a kind of hybrid optimization method, which combines direct optimization method and indirect optimization method. However, excessive series terms in OFS will lead to very time-consuming calculation. Thus, the computational efficiency is important considerations.
With regard to the sensitivity of the number of terms, we obtain a result that the selection of the number will influence the simulation obviously. Sometimes, we even cannot get the convergence solution. However, the shape-based method is proposed to get an approximate trajectory and provide an efficient initial guess for trajectory optimizers. So we only need to find a set of values which can satisfy the constraints.
6. Conclusion
A spacecraft’s initial trajectory design problem is studied in this paper. On the basis of the traditional Fourier series method, this paper proposes the optimized shape-based method using the finite Fourier series.
According to these simulation examples, we can get the following conclusions: Firstly, considering the first-order optimal necessary conditions, the optimized shape-based method can design the trajectory with minimal fuel consumption and offer initial guesses for detailed optimizers; secondly, different from other shape-based methods, the thrust direction of the spacecraft is not constrained to be tangential in the optimized shape-based method; thirdly, the shape-based method using finite Fourier series could be applied in the time-free case.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research is supported by the National Natural Science Foundation of China (Grant no. 11272255).