Lyapunovbased Lowthrust Optimal Orbit Transfer:
An approach in Cartesian coordinates
Abstract
This paper presents a simple approach to lowthrust optimalfuel and optimaltime transfer problems between two elliptic orbits using the Cartesian coordinates system. In this case, an orbit is described by its specific angular momentum and Laplace vectors with a free injection point. Trajectory optimization with the pseudospectral method and nonlinear programming are supported by the initial guess generated from the ChangChichkaMarsden Lyapunovbased transfer controller. This approach successfully solves several lowthrust optimal problems. Numerical results show that the Lyapunovbased initial guess overcomes the difficulty in optimization caused by the strong oscillation of variables in the Cartesian coordinates system. Furthermore, a comparison of the results shows that obtaining the optimal transfer solution through the polynomial approximation by utilizing Cartesian coordinates is easier than using orbital elements, which normally produce strongly nonlinear equations of motion. In this paper, the Earth’s oblateness and shadow effect are not taken into account.
keywords:
ChangChichkaMarsden Lyapunovbased transfer, Trajectory optimization, Cartesian coordinates1 Introduction
Three dimensional lowthrust optimal orbit transfer has attracted much inquiry focused on trajectory optimization using orbital elements. Due to strong nonlinearity of differential equations in Gaussian form with orbital elements, it is often difficult to obtain the optimal solution numerically in this system.
Some of the earliest works on the orbit transfer between neighboring elliptic orbits and on the transfer between coplanar and coaxial ellipses were presented by Edelbaum edelbaum1964optimum (); edelbaum1965optimum (). However, his elements as well as the Keplerian elements all contain singularity. To avoid the singularity, Arsenault arsenault1970 () firstly introduced the equinoctial elements. Broucke and Cefola broucke1972equinoctial () developed nonsingular equinoctial orbital elements using the Lagrange and Poisson brackets of Keplerian elements. Kechichian kechichian1996optimal () presented an application of these nonsingular elements to solve the minimum time rendezvous problem with constant acceleration. Chobotov chobotov2002orbital () considered more cases in minimum time transfer, including the comparison between exact solutions and approximate solutions obtained by the averaging technique. Gerffroy and Epenoy geffroy1997optimal () made further progress in both minimum time and fuel transfer problems using the averaging technique with the constraints of Earth’s oblateness and shadow effect taken into account. More recent works using the numerical averaging technique were presented by Tarzi and Speyer et al.tarzi2013fuel (), who provided a quick and accurate numerical approach for a wide range of transfers, including orbital perturbations such as Earth’s oblateness and shadow effect. Besides the strong nonlinearity, an additional disadvantage in using the equinoctial elements is the complexity of the equations in this coordinate system. With this system, Kepler’s equation must be solved by iteration to get the eccentric longitude at each integration step. Hence, equinoctial elements present challenges for trajectory optimization. Walker walker1985set (); walker1986set () put forth another important development in the study of equinoctial elements when he used the Stroboscopic method to modify orbital elements. He also altered differential equations into an approximative form, containing five dependent variables and one independent variable, as a means of achieving faster computation performance without solving Kepler’s equation. Roth roth1979higher () introduced the Stroboscopic method to obtain a higher order approximation for small perturbed dynamical systems, which depends on several slow variables and one fast variable. Haberkorn and Gergaud haberkorn2004low (); gergaud2006homotopy () investigated the application of the modified orbital elements by using the homotopy method. Cui et al. cui2008low () used sequential quadratic programming under their Lyapunov feedback control law, which was based on a function made up of modified elements, to obtain the optimalLyapunov solution without optimal transfer. Gao and Li gao2010optimization () made the similar work to optimize the Lyapunov function but never reached the optimal solution based on their Lyapunov control law.
An advanced technique using the Lyapunovbased controller to solve the lowthrust orbit transfer problem in Cartesian coordinates was introduced and rigorously proved by Chang et al. chang2002lyapunov (). This technique is based on the fact that a nondegenerate Keplerian orbit is uniquely described by its specific angular momentum and Laplace vectors. The resulting Lyapunov function provides an asymptotically stabilizing feedback controller, such that the target elliptic Keplerian orbit becomes a locally asymptotically stable periodic orbit. However, the Lyapunovbased transfer trajectory is not optimal in every sense. In this paper, the Lyapunovbased transfer presented in chang2002lyapunov () shall be called ChangChichkaMarsden (CCM)^{4}^{4}4Abbreviation for ChangChichkaMarsden transfer to distinguish it from any other Lyapunovbased transfers.
The motivation behind this paper is to use the CCM transfer trajectory as an initial guess for optimization in order to obtain the optimal transfer solution utilizing Cartesian coordinates. Using this method avoids the numerical disadvantages due to strong nonlinearity and complexity in the use of orbital elements.
This paper reviews the CCM transfer method chang2002lyapunov () in Section 2. In Section 3, a means to translate Keplerian elements into specific angular momentum and Laplace vectors is presented. Section 4 presents the proposed approach and the optimality (KKT)^{5}^{5}5Abbreviation for KarushKuhnTucker conditions for the minimum fuel consumption problem in Cartesian coordinates. Specifically, the ChebyshevGauss pseudospectral method is introduced to illustrate how the continuous optimal control problem can be reduced to a discretized nonlinear programming problem. Finally, in Section 5, numerical simulations are carried out to make detailed comparisons between the optimal results using Cartesian coordinates and those using orbital elements with the same initial guess. It shows that the use of Cartesian coordinates makes it easier to obtain the correct optimal solution.
2 ChangChichkaMarsden Transfer
This section summarizes and reviews the ChangChichkaMarsden (CCM) transfer in chang2002lyapunov (). This transfer employs Lyapunovbased controllers to achieve asymptotically stable transfers between elliptic orbits in a twobody problem.
2.1 TwoBody Problem
This paper assumes that the configuration space is . Let be the tangent space of , and be the coordinates on . The equations of motion are given by
(1) 
of which the solutions are regarded as the Keplerian orbits, where is the gravitational constant. The specific energy is defined by
Define by
where is the specific angular momentum vector and is the Laplace vector. The Laplace vector is related to the eccentricity vector as follows:
The three quantities , and satisfy the following two identities:
Define the sets
(2) 
The following Proposition is from chang2002lyapunov ().
Proposition 1.
The following hold:
1. is the union of all nondegenerate elliptic Keplerian orbits.
2. and .
3. The fiber consists of a unique (oriented) nondegenerate elliptic Keplerian orbit for each .
2.2 ChangChichkaMarsden Transfer
The equation of the motion with a specific control force is given by
(3) 
where is the control force. Define a metric on by
with an arbitrary parameter, and . Let be the open ball of radius centered at in the metric and its closure.
Let be the pair of the angular momentum and the eccentricity vector of a target elliptic orbit. Define a Lyapunov function on by
(4) 
Along the trajectory of (3) there is
Hence,
where
(5) 
Take a controller as follows:
(6) 
with an arbitrary function . Then,
(7) 
The following Proposition is proven in chang2002lyapunov () using LaSalle’s invariance principle (la1961stability, , pp.5859).
Proposition 2.
Let be the pair of the specific angular momentum and Laplace vectors of the target elliptic orbit. Take any closed ball of radius centered at contained in the set defined in (2). Then every trajectory starting in the subset of remains in that subset and asymptotically converges to the target elliptic orbit in the closedloop dynamics (3) with the control law in (6)
The choice of the parameter in the Lyapunov function plays a crucial role in determining the transfer trajectory chang2002lyapunov (). It determines the relative weighting between the two quadratic terms in the function in (4). With a small the shape of the trajectory will adjust to that of the target orbit first because a more weight is on . On the other hand, with a large , the normal direction of the trajectory plane will adjust to that of the target orbit plane first because a more weight is on . The parameter also determines the shape of the region of attraction since determines the shape of the ball in the metric . Additionally, the CCM transfer works well for parabolic transfer, although the success of the transfer is proven exclusively for elliptic orbits only in Proposition 2.
3 Transform of Orbital Elements
This section presents the transform of the six Keplerian elements to specific angular momentum and Laplace vectors for convenient reference. The state vector at periapsis to derive the transform in Cartesian coordinates with the Earth at the origin. Let
The periapsis of the orbit in the geocentric equatorial frame is determined by
In the perifocal frame,
The transformation matrix (curtis2005orbital, , p.174) from the perifocal frame into the geocentric equatorial frame is given by
The state vector in the geocentric equatorial frame is found by carrying out the matrix multiplications
Thus, the components of and are derived. Then using the identity , the specific angular momentum and Laplace vectors in the geocentric equatorial frame are computed as follows:
On equatorial orbits , they simplify to
On circular orbits , they simplify to
4 Optimal Orbit Transfer
The CCM transfer trajectory is used as an initial guess to support the trajectory optimization in the openloop system using the direct ChebyshevGauss pseudospectral transcription method and a nonlinear programming solver. The Cartesian coordinates and the modified orbital elements in optimization are compared.
4.1 Optimization in Cartesian Coordinates
Let denote the Cartesian coordinates in the geocentric equatorial frame. The minimum fuel consumption problem is given as follows:
with
where denotes the state vector, is the identity matrix and the control vector. The boundary conditions are given by
(8)  
(9)  
(10) 
For the minimum time problem, the cost function shall be used.
To reduce the continuous optimal control problem (OCP) into a discretized nonlinear programming (NLP) problem, the pseudospectral method is used with secondkind Chebyshev points.
The transformation to express the OCPs in the time interval is given by
Use Lagrange interpolation polynomials with points as follows:
where is the initial boundary point, and , , are the collocation points, which are the zeros of the secondkind Chebyshev polynomial as expressed below:
The weights of the ChebyshevGauss quadrature in this case are given by
Then the NLP problem can be obtained as (see (benson2005gauss, , pp.117–118), fahroo2002direct ())
Minimize  
Subject to  
where indicates the th collocating point. The augmented cost function with the constraints combined via Lagrange multipliers is given by
The remaining KarushKuhnTucker (KKT) conditions at the collocating points are given by
Now, the optimal control problem can be solved by using wellestablished NLP algorithms.
4.2 Modified Equinoctial Elements
The same notations as gergaud2006homotopy () are utilized here to describe the modified equinoctial orbit elements. The state variables are defined as
where the true longitude is the fast independent variable and the other five are slow dependent variables.
The control variables in RTN (RadialTangentialNormal) coordinates are defined with
Then
with
The system equation is given by
(11) 
The minimum fuel consumption orbit transfer problem can be written as
The optimality conditions expressed in terms of modified elements are not provided here since those presented above in Cartesian coordinates also apply to modified elements.
5 Numerical Results
Using an example, this section illustrates transfers from a lowEarth orbit (LEO) to a geosynchronous orbit (GEO) in terms of minimum transfer time and minimum fuel consumption. The numerical data used are from (chobotov2002orbital, , p.374), with the exception of the final time in the minimum fuel case. The initial point is given by
The target orbit is given by
with the control constraint km/sec and the initial time . In the minimum fuel case, the fixed final time is hr.
Canonical units are used in simulations, where 806.812 sec = 1 canonical time unit; 6378.140 km = 1 canonical distance unit; km/sec = 1 canonical acceleration unit; and the gravitational parameter . In the following, all units are canonical unless otherwise indicated. The initial and final conditions in Cartesian coordinates are given by
The initial and final conditions in modified elements are given by
with . In the minimum fuel case, the fixed final time is , but in the minimum time case the final time is free, so it must first be adjusted according to the final time guess of .
The CCM controller given in (6) is applied with the weighting and the function to obtain a transfer trajectory that provides an initial guess for optimization. Then TOMLAB/PROPT is utilized together with the pseudospectral method and SNOPT solver to optimize the trajectory on the MATLAB platform. To compare the different coordinate systems, all the collocating points are located in one phase. The optimal results are listed in Table 1, and the minimum fuel and minimum time transfer trajectories are shown in Fig. 1 and Fig. 2, respectively. For the minimum fuel case, Fig. 3 compares the force magnitude between the two coordinate systems presented above. Fig. 4 and Fig. 5 show the time history of , and , and Fig. 6 details the time history of the control force in the Cartesian coordinates system. In the time history plots, the dotted lines represent the initial guesses and the solid lines represent final trajectories.
Transfer  Switch  Collocating  Optimality  
Basis  Time (hr)  Times  (km/sec)  Points  Condition  
Method in this paper  
Cartesian  20.1703  6  4.9911  300  Satisfied  
20.1703  6  4.9938  350  Satisfied  
Min Fuel  20.1703  11  4.9945  300  Not Satisfied  
Modified  20.1703  12  5.1788  350  Not Satisfied  
20.1703  12  5.1772  400  Not Satisfied  
Min Time  Cartesian  16.1362  0  5.6873  300  Satisfied 
Modified  18.2205  0  6.4212  300  Satisfied  
Chobotov’s Result in (chobotov2002orbital, , pp.374376)  
Min Time  Equinoctial  16.2845  0  5.7451  *  * 
Narrower control constraints are then set by and to simulate the optimal fuel trajectory with the fixed final time and . The optimal results are listed in Table 2. The final transfer trajectory and the time histories of the state variables with are shown in Fig. 7 and Fig. 8.
Control  Transfer  Switch  Collocating  Optimality  

Basis  Constraints  Time (hr)  Times  (km/sec)  Points  Conditions  
20.1703  6  4.9938  350  
Min Fuel  Cartesian  40.3406  11  5.1185  400  Satisfied  
80.6812  31  5.1747  400 
In the minimum time case, the optimal result of transfer time obtained using the Cartesian coordinates system is similar to Chobotov’s in chobotov2002orbital (), while the optimized result obtained using modified elements is far from optimal. The mathematical description of the modified elements is an approximation from the equinoctial elements found by the stroboscopic method. Therefore, it fails to approximate the real optimal solution in the presented example. However, the equinoctial elements are not a good choice of a coordinate system for direct optimization since Kepler’s equation remains to be solved.
In the minimum fuel consumption case, the results from the use of Cartesian coordinates and those from the use of the modified elements with 300 collocating points in Table 1 are similar to each other. Even still, all the optimality conditions are not satisfied in the case where the modified elements are used. Fig. 3 shows six switching times of in the Cartesian coordinates system with 350 collocating points and eleven switching times of in the modified elements with 300 collocating points. This implies that using modified elements requires more collocating points and higher order polynomials to approximate the trajectory. Obtaining the optimal solution in this case increases the difficulty of computation noticeably. Even when modified elements are applied to find an optimal solution close to that which was obtained via Cartesian coordinates, this method would again fail to approximate the real optimal solution since the optimized control force sequences would be quite different from those in Cartesian coordinates. Moreover, problematic results can be identified with 350 and 400 collocating points from modified elements as illustrated in Table 1. However, by using Cartesian coordinates that contain weaker nonlinearity than orbital elements, obtaining the real optimal solution through polynomial approximation can be achieved with ease.
The time histories of in Fig. 5 and Fig. 8 indicate the importance of the CCM transfer as the initial guess of trajectory optimization in the Cartesian coordinates system. While the NLP solver is quite robust, without a supported Lyapnovbased initial guess, it would still be difficult for the solver to move those collocating points onto the optimal trajectory, due to the periodically changing sign and rapidly changing value of the state variables in the Cartesian coordinates system (conway2010spacecraft, , p.52).
Besides the strong nonlinearity in the equations of motion, the use of modified elements has several numerical disadvantages because its system, which is the righthand side of equation (11), is only an approximation. For a given lowthrust transfer trajectory, the time histories calculated from the righthand side of the differential equation (11) can not match the real time histories of the variables’ derivatives on lefthand side of (11). In particular, for stiff problems – as with long duration cases – where the trajectories are lacking high accuracy from the ODE solver, the error between the two sides of the equation will be large.
6 Conclusions
In this paper, a simple and effective approach has been proposed by employing the pseudospectral method, nonlinear programming, and the ChangChichkaMarsden transfer controller. Solutions to the freeinjection minimum fuel consumption and minimum time transfer problems have successfully been obtained by using Cartesian coordinates supported by a ChangChichkaMarsden transfer trajectory. The ChangChichkaMarsden transfer, as an initial guess, has revealed its usefulness for overcoming the numerical difficulties associated with the direct optimization caused by strong oscillation of state variables in the Cartesian coordinates system. The use of orbital elements increases the difficulty of optimization and fails to provide the optimal solution by the nonlinear programming and pseudospectral methods. However, by utilizing Cartesian coordinates instead of orbital elements, the optimal solutions are easy to obtain.
Two main advantages arise when Cartesian coordinates are used in the direct trajectory optimization. First, this option is the simplest way to describe the twobody system accurately without the optimal solutionrelated problems posed by singularity or approximation. Second, the weaker nonlinearity makes it easier to obtain the optimal solution, mainly due to the numerical approximation method used.
Future research on longer duration transfer described in Cartesian coordinates can be performed by developing a specific numerical method for quick and accurate computation with fewer collocating points. Additionally, the oblateness of the Earth and the shadow effect shall be taken into account.
7 Acknowledgments
This research was carried out in the Department of Applied Mathematics, University of Waterloo, and the School of Astronautics, Harbin Institute of Technology. The authors are indebted to Per Rutquist in TOMLAB for the advice on the use of the optimization software PROPT.
References

(1)
T. N. Edelbaum, ”Optimum LowThrust Rendezvous and Station Keeping,” Journal of Spacecraft and Rockets, Vol. 40, No. 6, 2003, pp. 960–965; Reprint from AIAA Journal, Vol. 2, No. 7, 1964, pp.1196–1201.
doi: 10.2514/2.7042. 
(2)
T. N. Edelbaum, ”Optimum PowerLimited Orbit Transfer in Strong Gravity Fields,” AIAA Journal, Vol. 3, No. 5, 1965, pp. 92–925.
doi:10.2514/3.3016. 
(3)
J. L. Arsenault, K. C. Ford, and P. E. Koskela, ”Orbit Determination using Analytic Partial Derivatives of Perturbed Motion,” AIAA Journal, Vol. 8, No. 1, 1970, pp. 4–12.
doi:10.2514/3.5597. 
(4)
R. A. Broucke, and P. J. Cefola, ”On the Equinoctial Orbit Elements,” Celestial Mechanics, Vol. 5, No. 3, 1972, pp. 303–310.
doi:10.1007/bf01228432. 
(5)
J. A. Kechichian, ”Optimal LowThrust Rendezvous using Equinoctial Orbit Elements,” Acta Astronautica, published online 18 Jun. 1999; Vol. 38, No. 1, 1996, pp. 1–14.
doi:10.1016/00945765(95)001212. 
(6)
V. A. Chobotov, Orbital Mechanics, 3rd ed., AIAA Education Series, AIAA, Inc., Reston, Virginia, 2002, pp. 372–389.
doi:10.2514/4.862250. 
(7)
S. Geffroy, and R. Epenoy, ”Optimal LowThrust Transfers with Constraints–Generalization of Aver aging Techniques,” Acta Astronautica, published online 17 Aug. 1998; Vol. 41, No. 3, 1997, pp. 133–149.
doi:10.1016/S00945765(97)002087. 
(8)
Z. Tarzi, J. L. Speyer, and R. E. Wirz, ”Fuel Optimum LowThrust Elliptic Transfer using Numerical Averaging,” Acta Astronautica, published online 16 Jan. 2013; Vol. 86, Jan. 2013, pp. 95–118.
doi:10.1016/j.actaastro.2013.01.003. 
(9)
M. J. H. Walker, B. Ireland, and J. Owens, ”A Set Modified Equinoctial Orbit Elements,” Celestial Mechanics, Vol. 36, No. 4, 1985, pp. 409–419.
doi:10.1007/bf01227493. 
(10)
M. J. H. Walker, ”A Set of Modified Equinoctial Orbit Elements,” Celestial Mechanics, Vol. 38, No. 4, 1986, pp. 391–392.
doi:10.1007/bf01238929. 
(11)
E. A. Roth, ”On the HigherOrder Stroboscopic Method,” Zeitschrift für angewandte Mathematik und Physik ZAMP, Vol. 30, No. 2, 1979, pp. 315–325.
doi:10.1007/bf01601943. 
(12)
T. Haberkorn, P. Martinon, and J. Gergaud, ”Low Thrust MinimumFuel Orbital Transfer: A Homotopic Approach,” Journal of Guidance, Control, and Dynamics, Vol. 27, No. 6, 2004, pp. 1046–1060.
doi:10.2514/1.4022. 
(13)
J. Gergaud, and T. Haberkorn, ”Homotopy Method for Minimum Consumption Orbit Transfer Problem,” ESAIM: Control, Optimisation and Calculus of Variations, published online 22 Feb. 2011; Vol. 12, No. 2, 2006, pp. 294–310.
doi:10.1051/cocv:2006003. 
(14)
P. Cui, Y. Ren, and E. Luan, ”LowThrust, MultiRevolution Orbit Transfer under the Constraint of a Switch Function without Prior Information,” Transactions of the Japan Society for Aeronautical and Space Sciences, published online 20 Feb. 2008; Vol. 50, No. 170, 2008, pp. 240–245.
doi:10.2322/tjsass.50.240. 
(15)
Y. Gao, and X. Li, ”Optimization of Lowthrust Manyrevolution Transfers and Lyapunovbased Guidance”, Acta Astronautica, Vol. 66, No. 1–2, 2010, pp. 117–129.
doi:10.1016/j.actaastro.2009.05.013. 
(16)
D. E. Chang, D. F. Chichka, and J. E. Marsden, ”LyapunovBased Transfer Between Elliptic Keplerian Orbits,” Discrete and Continuous Dynamical Systems Series B, published online Nov. 2001; Vol. 2, No. 1, 2002, pp. 57–67.
doi:10.3934/dcdsb.2002.2.57.  (17) J. P. La Salle, and S. Lefschetz, Stability by Liapunov’s Direct Method with Applications, Vol. 4, Mathematics in Science and Engineering Series, Academic Press, New York, 1961, pp. 58–59.
 (18) H. D. Curtis, Orbital Mechanics for Engineering Students, 1st ed., Elsevier Aerospace Engineering Series, Elsevier ButterworthHeinemann, Oxford, 2005, p.174.
 (19) D. Benson, A Gauss Pseudospectral Transcription for Optimal Control, Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA, Feb. 2005.

(20)
F. Fahroo, and I. M. Ross, ”Direct Trajectory Optimization by a Chebyshev Pseudospectral Method,” Journal of Guidance, Control, and Dynamics, Vol. 25, No. 1, 2002, pp. 160–166.
doi:10.2514/2.4862.  (21) B. A. Conway, Spacecraft Trajectory Optimization, Cambridge Aerospace Series, Cambridge University Press, New York, 2010, p. 52.
Hantian Zhang is a current undergraduate student at Department of Aerospace Engineering and Mechanics, School of Astronautics, Harbin Institute of Technology. He expected to receive his B.Eng. from Harbin Institute of Technology in 2014. He worked as a research assistant with Prof. Dong Eui Chang at Department of Applied Mathematics in University of Waterloo from January to March in 2013. His research interests include orbital dynamics and control, trajectory optimization, nonsmooth mechanics, and dynamical systems.  
Dong Eui Chang received his B.S. and M.S. from Seoul National University in 1994 and 1997, and his Ph.D. in Control & Dynamical Systems (CDS) from the California Institute of Technology (Caltech). He is currently an Associate Professor in the Department of Applied Mathematics at the University of Waterloo. His research interests lie in control, mechanics and various engineering applications.  
Qingjie Cao received his B.S. and M.S. in Mathematics from Qufu Normal University in 1982 and 1985, and his Ph.D. in Mechanical Engineering from Tianjin University in 1993. Before he joined Harbin Institute of Technology, he served as a Professor at School of Mathematics in Shandong University, and Research Fellows in University of Aberdeen and University of Liverpool. Currently, he is a Professor at School of Astronautics in Harbin Institute of Technology, and the director of the Centre for Nonlinear Dynamics Research. He introduced the SD oscillator and SD attractor in mechanical systems. 