Research Article  Open Access
A ShapeBased Method for Continuous LowThrust Trajectory Design between Circular Coplanar Orbits
Abstract
The shapebased 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 shapebased method using the finite Fourier series is proposed in this paper. Taking the flight timefixed case and the timefree case into account, respectively, the optimized shapebased method, which considers the firstorder 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 shapebased method using the finite Fourier series, in which the thrust direction is constrained to be tangential, is developed for the timefree case in this paper. The EarthMars case and the LEOGEO case are used to verify the optimized shapebased method’s feasibility for timefixed and timefree continuous lowthrust trajectory design between circular coplanar orbits, respectively. The optimized shapedbased method can design a lower cost trajectory.
1. Introduction
Recently, continuous lowthrust trajectory design and optimization are becoming increasingly popular [1, 2], although they are very challenging and timeconsuming. What is particular is that the continuous lowthrust 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 shapebased (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 twodimensional (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 lowthrust trajectory design procedure. Cui et al. [7] proposed a new search approach algorithm for the launch window of lowthrust gravityassist 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 spiralbased (LS) nonKeplerian orbit. The feasibility and essential characteristics of the LS nonKeplerian 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 abovementioned disadvantages, Wall and Conway [9, 10] developed a 2D inverse polynomial (IP) method. The fifthorder IP method can be used to design the transfer trajectory for the timefree case, while the sixthorder IP method is designed for the timefixed 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 lowthrust transfer trajectory. Considering thrust and radius constraints, Wang et al. [12] proposed a modified IP method for both the timefree transfer case and the timefixed 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 lowthrust trajectory design between elliptical orbits, Xie et al. [13] constructed a new shape function about two semimajoraxis parameters, which are polynomials in the polar angle. With tangential thrust, a fifthorder and sixthorder method is designed for timefree and timefixed cases as in the IP method.
Taheri and Abdelkhalik [14] and Abdelkhalik and Taheri [15] proposed a new shapebased 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 lowthrust trajectory design method called hodographicshaping (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 firstorder 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 shapebased method using finite Fourier series, which can easily overcome the abovementioned shortcomings, is proposed in this paper. Regarding spacecraft threedimensional (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 timefixed rendezvous case and timefree transfer case, respectively, the firstorder 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 shapebased method’s feasibility for continuous lowthrust 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 ShapeBased Method Using Fourier Series
3.1. Performance Index and Boundary Conditions
Usually the performance index for the lowthrust 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 firstorder or secondorder derivative of a certain variable, is more reasonable.
3.2. The Flight TimeFixed Case
3.2.1. Traditional Fourier Series (TFS) Method
In the timefixed case, it is assumed that the thrust is aligned along or against the velocity vector; that is, , where . The TFS method for the timefixed 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 timefixed 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 timefixed with no limitation on thrust direction, (4), (5) and (22) become an NLP problem about Fourier coefficients.
3.3. The Flight TimeFree Case
3.3.1. TFS Method
In the flight timefree 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 timefree 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 timefree 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 timefree 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 lowthrust 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 timefixed and timefree cases, the continuous lowthrust EarthMars rendezvous and LEOGEO 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 EarthMars Rendezvous
The FS method for the timefixed case can be applied to the interplanetary exploration, asteroid deflection, rendezvous and docking missions, and so on. In the timefixed case, we design the continuous lowthrust EarthMars 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 EarthMars 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 lowthrust 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 LEOGEO Transfer
The FS method for the timefree case can also be applied for a spacecraft orbit raise and maneuver mission. For the timefree case, the LEOGEO 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 timefree case, and their firstorder derivatives with respect to polar angle shall be calculated as new boundary conditions, as shown in
Figure 8 gives the spacecraft’s LEOGEO 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 LEOGEO case.
Figures 8 and 10 prove that the shapebased trajectory design method using the finite Fourier series can be used in the flight timefree case. This paper extends the scope of application of the Fourier series method to the flight timefree case. Figure 10 shows that, in the flight timefree 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 EarthMars transfer in the flight timefixed 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 EarthMars 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 timefixed case, the OFS method in the timefree 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 timeconsuming 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 shapebased 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 shapebased method using the finite Fourier series.
According to these simulation examples, we can get the following conclusions: Firstly, considering the firstorder optimal necessary conditions, the optimized shapebased method can design the trajectory with minimal fuel consumption and offer initial guesses for detailed optimizers; secondly, different from other shapebased methods, the thrust direction of the spacecraft is not constrained to be tangential in the optimized shapebased method; thirdly, the shapebased method using finite Fourier series could be applied in the timefree 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).
References
 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  Zentralblatt MATH
 R. J. McKay, M. MacDonald, J. Biggs, and C. McInnes, “Survey of highlynonKeplerian orbits with lowthrust propulsion,” Journal of Guidance, Control, and Dynamics, vol. 34, no. 3, pp. 645–666, 2011. View at: Publisher Site  Google Scholar
 S. Li, Y. Zhu, and Y. Wang, “Rapid design and optimization of lowthrust rendezvous/interception trajectory for asteroid deflection missions,” Advances in Space Research, vol. 53, no. 4, pp. 696–707, 2014. View at: Publisher Site  Google Scholar
 A. E. Petropoulos, Shapebased algorithm for the automated design of lowthrust, gravity assist trajectories [Ph.D. thesis], Purdue University, West Lafayette, Ind, USA, 2001.
 A. E. Petropoulos and J. M. Longuski, “Shapebased algorithm for automated design of lowthrust, gravityassist trajectories,” Journal of Spacecraft and Rockets, vol. 41, no. 5, pp. 787–796, 2004. View at: Publisher Site  Google Scholar
 D. Izzo, “Lambert's problem for exponential sinusoids,” Journal of Guidance, Control, and Dynamics, vol. 29, no. 5, pp. 1242–1245, 2006. View at: Publisher Site  Google Scholar
 P.Y. Cui, H.B. Shang, and E.J. Luan, “A fast search algorithm for launch window of interplanetary lowthrust exploration mission,” Journal of Astronautics, vol. 29, no. 1, pp. 40–45, 2008 (Chinese). View at: Google Scholar
 L.L. Zheng, J.P. Yuan, and Z.X. Zhu, “Logarithmic spiralbased nonKeplerian orbit design,” Journal of Astronautics, vol. 31, no. 9, pp. 2075–2081, 2010 (Chinese). View at: Publisher Site  Google Scholar
 B. J. Wall, “Shapebased approximation method for lowthrust trajectory optimization,” in Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, Hawaii, USA, August 2008. View at: 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–102, 2009. View at: Publisher Site  Google Scholar
 H. Shang, P. Cui, D. Qiao, and R. Xu, “Lambert solution and application for interplanetary lowthrust trajectories,” Acta Aeronautica et Astronautica Sinica, vol. 31, no. 9, pp. 1752–1757, 2010 (Chinese). View at: Google Scholar
 D. Wang, G. Zhang, and X. Cao, “Modified inversepolynomial shaping approach with thrust and radius constraints,” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 229, no. 13, pp. 2506–2518, 2015. View at: Publisher Site  Google Scholar
 C. Xie, G. Zhang, and Y. Zhang, “Simple shaping approximation for lowthrust trajectories between coplanar elliptical orbits,” Journal of Guidance, Control, and Dynamics, vol. 38, no. 12, pp. 2448–2455, 2015. View at: Publisher Site  Google Scholar
 E. Taheri and O. Abdelkhalik, “Shapebased 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
 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
 D. J. Gondelach, R. Noomen, and D. Geller, “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
 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, Rapid space trajectory generation using a Fourier series shapebased approach [Ph.D. thesis], Michigan Technological University, 2014.
 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
Copyright
Copyright © 2017 Qun Fang 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.