Lyapunov-based Low-thrust Optimal Orbit Transfer:An approach in Cartesian coordinates

Lyapunov-based Low-thrust Optimal Orbit Transfer:
An approach in Cartesian coordinates

Hantian Zhang111Undergraduate Student; School of Astronautics, Harbin Institute of Technology, Harbin, China, 150001 Dong Eui Chang222Associate Professor; Department of Applied Mathematics, University of Waterloo, Waterloo, ON, Canada, N2L 3G1 Qingjie Cao333Professor;

This paper presents a simple approach to low-thrust optimal-fuel and optimal-time 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 Chang-Chichka-Marsden Lyapunov-based transfer controller. This approach successfully solves several low-thrust optimal problems. Numerical results show that the Lyapunov-based 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.

Chang-Chichka-Marsden Lyapunov-based transfer, Trajectory optimization, Cartesian coordinates

1 Introduction

Three dimensional low-thrust 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 optimal-Lyapunov 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 Lyapunov-based controller to solve the low-thrust 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 non-degenerate 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 Lyapunov-based transfer trajectory is not optimal in every sense. In this paper, the Lyapunov-based transfer presented in chang2002lyapunov () shall be called Chang-Chichka-Marsden (CCM)444Abbreviation for Chang-Chichka-Marsden transfer to distinguish it from any other Lyapunov-based 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)555Abbreviation for Karush-Kuhn-Tucker conditions for the minimum fuel consumption problem in Cartesian coordinates. Specifically, the Chebyshev-Gauss 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 Chang-Chichka-Marsden Transfer

This section summarizes and reviews the Chang-Chichka-Marsden (CCM) transfer in chang2002lyapunov (). This transfer employs Lyapunov-based controllers to achieve asymptotically stable transfers between elliptic orbits in a two-body problem.

2.1 Two-Body 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


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


The following Proposition is from chang2002lyapunov ().

Proposition 1.

The following hold:
1. is the union of all non-degenerate elliptic Keplerian orbits.
2. and .
3. The fiber consists of a unique (oriented) non-degenerate elliptic Keplerian orbit for each .

2.2 Chang-Chichka-Marsden Transfer

The equation of the motion with a specific control force is given by


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


Along the trajectory of (3) there is




Take a controller as follows:


with an arbitrary function . Then,


The following Proposition is proven in chang2002lyapunov () using LaSalle’s invariance principle (la1961stability, , pp.58-59).

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 closed-loop 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 open-loop system using the direct Chebyshev-Gauss 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:


where denotes the state vector, is the identity matrix and the control vector. The boundary conditions are given by


For the minimum time problem, the cost function shall be used.

To reduce the continuous optimal control problem (OCP) into a discretized non-linear programming (NLP) problem, the pseudospectral method is used with second-kind 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 second-kind Chebyshev polynomial as expressed below:

The weights of the Chebyshev-Gauss quadrature in this case are given by

Then the NLP problem can be obtained as (see (benson2005gauss, , pp.117–118), fahroo2002direct ())

Subject to

where indicates the th collocating point. The augmented cost function with the constraints combined via Lagrange multipliers is given by

The remaining Karush-Kuhn-Tucker (KKT) conditions at the collocating points are given by

Now, the optimal control problem can be solved by using well-established 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 (Radial-Tangential-Normal) coordinates are defined with



The system equation is given by


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 low-Earth 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.374-376)
Min Time Equinoctial 16.2845 0 5.7451 * *
Table 1: Comparisons of optimal results in different basis
Figure 1: Minimum fuel transfer (Cartesian)
Figure 2: Minimum time transfer (Cartesian)
(a) (350 points, Cartesian)
(b) (300 points, Modified)
Figure 3: Comparison of control force magnitude time history of min fuel
Figure 4: time history of min fuel (Cartesian)
Figure 5: time history of min fuel (Cartesian)
Figure 6: Control force time history of min fuel (Cartesian)

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
Table 2: Optimal results of low-thrust minimum fuel problem
Figure 7: Minimum fuel transfer (Cartesian, )
Figure 8: time history of min fuel (Cartesian, )

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 Lyapnov-based 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 right-hand side of equation (11), is only an approximation. For a given low-thrust transfer trajectory, the time histories calculated from the right-hand side of the differential equation (11) can not match the real time histories of the variables’ derivatives on left-hand 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 Chang-Chichka-Marsden transfer controller. Solutions to the free-injection minimum fuel consumption and minimum time transfer problems have successfully been obtained by using Cartesian coordinates supported by a Chang-Chichka-Marsden transfer trajectory. The Chang-Chichka-Marsden 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 two-body system accurately without the optimal solution-related 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.


  • (1) T. N. Edelbaum, ”Optimum Low-Thrust 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 Power-Limited Orbit Transfer in Strong Gravity Fields,” AIAA Journal, Vol. 3, No. 5, 1965, pp. 92–925.
  • (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.
  • (4) R. A. Broucke, and P. J. Cefola, ”On the Equinoctial Orbit Elements,” Celestial Mechanics, Vol. 5, No. 3, 1972, pp. 303–310.
  • (5) J. A. Kechichian, ”Optimal Low-Thrust Rendezvous using Equinoctial Orbit Elements,” Acta Astronautica, published online 18 Jun. 1999; Vol. 38, No. 1, 1996, pp. 1–14.
  • (6) V. A. Chobotov, Orbital Mechanics, 3rd ed., AIAA Education Series, AIAA, Inc., Reston, Virginia, 2002, pp. 372–389.
  • (7) S. Geffroy, and R. Epenoy, ”Optimal Low-Thrust Transfers with Constraints–Generalization of Aver- aging Techniques,” Acta Astronautica, published online 17 Aug. 1998; Vol. 41, No. 3, 1997, pp. 133–149.
  • (8) Z. Tarzi, J. L. Speyer, and R. E. Wirz, ”Fuel Optimum Low-Thrust Elliptic Transfer using Numerical Averaging,” Acta Astronautica, published online 16 Jan. 2013; Vol. 86, Jan. 2013, pp. 95–118.
  • (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.
  • (10) M. J. H. Walker, ”A Set of Modified Equinoctial Orbit Elements,” Celestial Mechanics, Vol. 38, No. 4, 1986, pp. 391–392.
  • (11) E. A. Roth, ”On the Higher-Order Stroboscopic Method,” Zeitschrift für angewandte Mathematik und Physik ZAMP, Vol. 30, No. 2, 1979, pp. 315–325.
  • (12) T. Haberkorn, P. Martinon, and J. Gergaud, ”Low Thrust Minimum-Fuel Orbital Transfer: A Homotopic Approach,” Journal of Guidance, Control, and Dynamics, Vol. 27, No. 6, 2004, pp. 1046–1060.
  • (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.
  • (14) P. Cui, Y. Ren, and E. Luan, ”Low-Thrust, Multi-Revolution 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.
  • (15) Y. Gao, and X. Li, ”Optimization of Low-thrust Many-revolution Transfers and Lyapunov-based Guidance”, Acta Astronautica, Vol. 66, No. 1–2, 2010, pp. 117–129.
  • (16) D. E. Chang, D. F. Chichka, and J. E. Marsden, ”Lyapunov-Based Transfer Between Elliptic Keplerian Orbits,” Discrete and Continuous Dynamical Systems Series B, published online Nov. 2001; Vol. 2, No. 1, 2002, pp. 57–67.
  • (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 Butterworth-Heinemann, 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.
  • (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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description