Available online at www.sciencedirect.com

ScienceDirect

Chinese Journal of Aeronautics 22(2009) 56-69

Chinese Journal of Aeronautics

www.elsevier.com/locate/cj a

Systematic Direct Approach for Optimizing Continuous-thrust

Earth-orbit Transfers

Gao Yanga'*, Li Weiqib

"Academy of Opto-Electronics, Chinese Academy of Sciences, Beijing 100080, China

bSchool of Automation Science and Electrical Engineering, Beijing University of Aeronautics and Astronautics, Beijing 100191, China

Received 21 March 2008; accepted 19 May 2008

Abstract

This article presents a systematic direct approach to carry out effective optimization of a wide range of continuous-thrust Earth-orbit transfers with intermediate-level thrust acceleration, including minimum-time (with a single burn arc) and minimum-fuel (with multiple burn arcs) transfers. With direct control parameterization, in which the control steering programs of burn arcs are interpolated through a finite number of nodes, the optimal control problem is converted into the parameter optimization problem to be solved by nonlinear programming. An inertial coordinate transformation strategy is proposed to produce general simple terminal constraints for transfers to inclined and eccentric orbits. This strategy eliminates possible singularities when the spacecraft transfers to the retrograde equatorial orbit and also results in better convergence of optimization iterations. Multiple shooting is used in order to further improve the convergence robustness, and strategies to allocate multiple-shooting variables for certain types of transfers are presented. Numerical examples include single-burn minimum-time transfers, multiple-burn minimum-fuel transfers with up to 12 burn arcs as well as fixed-time transfers with variable specific impulses. All the converged solutions are obtained from simply defined initial guesses. The proposed direct approach is briefly compared with other methods. Finally, we discussed the potential onboard guidance scheme using model predictive control.

Keywords: low-thrust; orbital transfer; optimization; multiple shooting; variable specific impulse; guidance scheme; model predictive control

1. Introduction

For decades, many researchers have been studying optimal orbit transfers by using continuous-thrust or low-thrust propulsion (usually referred to as electric propulsion), which is viewed as a prospective alternative for near-term space missions. The spacecraft, Deep Space 1, witnessed the first successful use of electric propulsion for an interplanetary mission[1]. Recently, the practical operation of solar electric propulsion was also investigated through the Smart-1 lunar mission[2]. Continuous-thrust propulsion is usually related to high specific impulse values, resulting in prolonged engine operation period during transfers.

Corresponding author. Tel.: +86-10-62582810.

E-mail address: gaoy@aoe.ac.cn

Foundation items: National Natural Science Foundation of China (10603005); Foundation of President of the Academy of Opto-Electro-nics ( A0E-CX-200601)

1000-9361/$ - see front matter © 2009 Elsevier Ltd. All rights reserved. doi: 10.1016/S1000-9361(08)60069-2

Orbit transfers by using continuous-thrust propulsion require continuous control profiles that cannot be approximated by impulsive velocity changes, which definitely poses a challenge to the task of designing transfer trajectories with optimal performances. The continuous-thrust electric propulsion is specifically referred to as low-thrust propulsion in the most research documents. However, the current technology of electric propulsion shows that the low-thrust propulsion can generate thrust as high as chemical propulsion. Therefore, we use the term "continuous-thrust" in this study to replace the term "low-thrust".

Up to date, researchers have been mainly using numerical methods such as indirect, direct, and so-called hybrid methods in continuous-thrust trajectory optimization. J. A. Kechichian[3] solved optimal Earth-orbit intermediate acceleration transfers through formulation of the two-point boundary-value problem (TPBVP). This is usually termed as indirect method characterized by the fact that the solution of the TPBVP is highly sensitive to the initial guesses of the

costate variables. This sensitivity can be alleviated if we resort to direct methods. A. L. Herman, et al.[4] developed a collocation method on the base of high-order Gauss-Lobatto quadratic rules, which had been used by A. L. Herman and D. B. Spencer[5] to deal with a wide range of optimal Earth-orbit transfers by using nonlinear programming (NLP). The high-order curve fitting rule could reduce the dimension of NLP problems in comparison with the traditional HermiteSimpson method. W. A. Scheel, et al.[6] introduced a parallel Runge-Kutta method to solve coplanar continuous-thrust transfer problems, leading to reduced number of optimization variables compared with the collocation method. There is a kind of hybrid methods that use costate equations to govern the optimal control without considering constraints of the switching function and transversality condition derived from the TPBVP. With the hybrid method, K. P. Zondervan, et al.[7] obtained a 3-burn transfer with a large plane change; Y. Gao, et al.[8] effectively solved a wide range of interplanetary transfer problems; M. R. Ilgen[9] obtained continuous-thrust Earth-orbit transfers. For very-low-continuous-thrust many-revolution transfers, the orbital averaging technique is widely used[10-12]. However, this sort of transfers is not considered in this article. Transfers with multiple-burn arcs at the expense of prolonged transfer time were investigated in Refs. [13-16]. The multiple-burn scheme is beneficial for fuel-saving, especially for transfers with a small number of revolutions. Recently, though nongradient-based methods[17-18] have found some application in trajectory optimization, their capability and efficiency are still limited in comparison with gradient-based methods.

This article presents a systematic direct approach that uses direct control parameterization to solve a wide range of Earth-orbit transfers with thrust-to-weight ratios on the order above 103. This sort of orbit transfers results in a small number of transfer revolutions and excludes the use of orbital averaging. The proposed approach was classified as a direct-multiple-shooting method by D. G. Hull[19] to solve general optimal control problems. In this study, a coordinate transformation strategy is introduced to produce general simple expressions of terminal constraints beneficial to improve the convergence robustness and also remove possible singularities of equinoctial elements in the retrograde equatorial orbit. The Earth J2 perturbations in the transformed inertial coordinate are derived and included in trajectory optimization. Furthermore, multiple shooting, categorized as the multiple shooting for thrusting arcs and multiple shooting for burn/coast sequences, is utilized to further improve the convergence robustness of NLP problems. How to allocate multiple-shooting variables is described for different types of Earth-orbit transfer problems. The systematic direct approach is able to effectively solve optimal single-burn minimum-time and multiple-burn minimum-fuel transfer problems by using the gradi-

ent-based optimization method. In addition, the capability of changeable specific impulses can be well accommodated and the effects of modulating specific impulses for Earth-orbit transfers are presented. In terms of excellent convergence of optimization problems, we finally discussed the potential guidance scheme using model predictive control.

2. System Dynamics and Direct Control Parameterization

A set of equinoctial elements introduced in Refs. [20-21] is utilized to eliminate singularities when either eccentricity or inclination equals zero. The equinoctial elements (p, f g, h, k, L) can be obtained in terms of six classical orbital elements: p = a(1 - e2), f = ecos(®+^) , g = esin(®+^) , h = tan(i/2)« cos C2 , k = tan(i/2)sin f2 , and L = Q + m+Q ,

where a is semimajor axis, e eccentricity, i inclination, Q right ascension of ascending node (RAAN), a> argument of perigee, and 6 true anomaly. Note that the equinoctial elements accommodate to all possible two-body conic orbits except in the case where i = 180°. The equations of motion of a spacecraft expressed by the equinoctial elements are

P - ^fe (1)

f = P i fr sin L + [(w +1) cos L + f ]-fe-

V<" I w

(h sin L - k cos L) gfl (2)

g = Jp fr cos L + [(w + 1)sin L + g]—fe +

(hsin L -kcos L)ffl (3)

h piUjL cos L (4)

\ H 2w

k = psf sin L (5)

\ ju 2w

L = J^pf —I +—. P(h sin L - kcos L) fn (6) IP) wV M

m = -T /(g0/,p) (7)

where w = 1 + f cos L + g sin L , s2 = 1 + h2 + k2, T is the thrust magnitude, m the spacecraft mass, g0 (= 9.806 65 m/s2) the Earth sea-level gravitational acceleration, and Ip the engine specific impulse, and H (=398 601 m /s2) the Earth gravitational parameter.

Three acceleration components fr, fe, and fn in the local vertical and local horizon (LVLH) frame are de-

fined by (the subscripts r, 0, and n represents the radial, circumferential, and normal directions, respectively)

[ fr fe fn ]T = - « + f

where /p=[fj2,r fj2fi fj2,n] is the perturbation acceleration with only J2 perturbations included in this work, and a is the direction unit vector of thrust acceleration, which can be expressed in terms of pitch and yaw steering angles in the LVLH frame:

a = [ar ae an ]T = [sin £cos ^ cos 8cos$ sin (9)

where the pitch angle 8 is measured from the local horizon to the projection of the thrust vector onto the orbital plane, and the yaw angle ^ from the orbital

plane to the thrust vector. The thrust magnitude is modeled by

T = (T / W)0 m0g0 (10)

where (T / W)0 is the initial thrust-to-weight ratio, representing a single physical variable. The initial spacecraft mass is denoted by m0.

The time histories of control steering are represented by a finite number of nodes, and the continuous control steering is obtained by using linear or cubic spline interpolation through these nodes. The optimal control problem is converted into the parameter optimization one which is in turn solved by NLP—se-quential quadratic programming (SQP)[22]. This method, termed direct control parameterization, was once used by C. A. Kluever[23] to solve a wide range of interplanetary orbit transfers. In this study, the direction

cosines of control direction or steering angles are parameterized into discrete nodes. The typical initial guess of the control steering nodes is defined by

a = [ar ae an]T=[0 1 0]T or 8 = 0, <p = 0 (11)

Eq.(11) yields a circumferential control steering as the initial guess for trajectory optimization.

3. Inertial Coordinate Transformation Strategy

The J2000 Earth equatorial frame was defined by the mean orientation of the Earth's equator and ecliptic orbit at the beginning of the year 2000. The Oxy plane is the plane of the mean Earth's equator. The line formed by intersection of the ecliptic orbit plane and the Earth's equatorial plane defines the x axis. The line starting from the sun to the center of the Earth on the first day of autumn defines the positive direction of the x axis called the vernal equinox. The z axis perpendicular to the Oxy plane points to the north. The y axis constitutes Cartesian coordinate according to the right-handed principle. Assuming that the original J2000 Earth equatorial inertial coordinate is denoted by Oxyz, an arbitrary inertial coordinate could be obtained by performing the following transformations in terms of three consecutive Euler-angle rotations: ® to rotate Oxyz by an angle a about the z axis to obtain Ox1y1z1; ® to rotate Ox1y1z1 by an angle p about the x1 axis to obtain Ox2y2z2; ® to rotate Ox2y2z2 by an angle y about the z2 axis to obtain Oxyz . This series of rotations is denoted by a rotation matrix A written as follows

cos ycos a-cos ^sin ysin a cos ysin a + sin ^cos ^cos a sin ysin ß - sin ycos a-cos ^cos ^sin a - sin ^sin a+ cos ^cos ^cos a cos ^sin ß

sin ^sin a

- sin ^cos a

The transformation from the coordinates Oxyz to Oxyz is

x x x x x

y = A y > y = A-1 y AT y

z z z S z

In the transformed inertial coordinate (denoted by Oxyz), the equations of motion (Eqs.(1)-(7)) remain

unchanged, but the six equinoctial elements are all referenced to Oxyz not to Oxyz (the variables with are meant being referenced in Oxyz). If the classical orbital elements of the target orbit are af, ef, if, i2t, ai, and dt in Oxyz, we can obtain the terminal constraints of equinoctial elements expressed in Oxyz with three Euler-angle rotations: a=Q{, p=if, and y = a{.

pf = af(1 -ef2),ff = ef,gf = 0,hf = k{ = 0,Lf =0f (14)

In Oxyz , the x axis points to the direction of the eccentricity vector; the orbit plane of the target orbit (defined in Oxyz) is the Oxy! plane; the z axis is perpendicular to the Oxy! plane. If the target orbit is a circular orbit (ef = 0 ), the first two rotations are enough without considering the third because of the undefined argument of perigee. If the target orbit is a circular equatorial orbit (ef = 0 and if = 0), no transformation is needed. Since the equations of motion are integrated over the equinoctial elements, experience about solving optimal transfer problems with SQP[22] shows that the constraints expressed in Eq.(14) are more robust for convergence than expressed by the classical orbital elements that are transformed non-linearly from the terminal equinoctial elements. Another advantage lies in avoiding singularity (i = 180°) by choosing appropriate values of the rotation angles for orbit transfers to the retrograde equatorial orbit,

which will be shown in Section 6 without re-defining the equinoctial elements.

Before starting optimization, the initial orbit elements (a0, e0, i0, a>0, and 00) defined in Oxyz

are transformed into a0, e0, i0, /30, a>0 , and 90 (with a=n{, p=i{, and y = a{) in Oxyz where the trajectory optimization is performed. In fact, a0, e0,

and 90 remain the same as a0, e0

i0, £20, ®0 might be different from i0, /20, a>0. After achieving optimal solutions in , an inverse

transformation is needed to obtain the states referenced in Oxyz. Note that the pitch and yaw angles (see Eq. (9)) expressed in the LVLH frame (either in Oxyz or in Oxyz ) are not affected by the inertial coordinate transformation.

4. J2 Perturbations Expressed in the Transformed Coordinate

The potential function of J2 perturbations defined in Oxyz can be written into

U = 4 R2 J21(3 sin2 0 -1), sin 0 = - (15) 2 r 2 r

where Re is the radius of Earth, and J2 the second-order zonal harmonic coefficient. In the new coordinate Oxyz, we can obtain Eq.(16) according to Eqs.(12)-(13).

z = alxc + a2y + a3z , r = r (16)

where al = sin ^sin p, a2 = cos ^sin p, a3 = cos p,

and r = *Jx2 + y2 + z2 . With sin ^ = sin i sin m , where m = rn + d , denoting the longitude angle, the potential function J2 perturbations, expressed by the classical orbit elements in Oxyz , becomes

Uj =4 R J 2 •—13[ aj(cos/2 cos sr -2 r 21

sin /2 sin sr cos i) + a2(sin /2 cos sr +

cos /2sinsrcos i) + a3sinsrcos i]2 -1} (17)

The acceleration components of the J2 perturbations in the LVLH frame can be obtained as follows

J = Re J 2(3 Sin2 ¿-1)

1 U n 2

f,« = "4Re J2*

"2'~ r dm r' 6 sin ^[aj (- cos /2 sin sr - sin Q cos sr cos i) -a2 (- sin /2 sin sj + cos /2 cos sr cos i) +

a3cossr sin i)] (19)

1 3U, u ,

•/J,,n _ ~ • ~ -v _ _4'ne,y2

2 r sin TU 01 r

6 sin ^(aj sin /2 sin i - a2 sin /2 sin i + a3 cos i) (20) where

sin tj) = al (cos /2 cos sr - sin /2 sin sr cos i) + a2 (sin /2 cos sr + cos /2 sin sr cos i) +

a3sinsr sin i (21)

The equinoctial elements can be used to express the acceleration of the J2 effects based on Eqs.(22)-(24).

cos n = h /V h2 + k2

sin ^ = k/Vh2 + k2 }

cos sr = (h cos L + k sin L)/V h2 + k2

sin ^ = (h sin L - k cos L) / Vh2 + k2 J

cosi = (1 - h2 - k 2)/(1 + h2 + k2) I

sin i = 2>/h2 + k2/(1 + h2 + k2) J Here we can readily obtain

1 —(

r' 2 sin i

J = Re2J2 • 1(3 sin2 1)

,, 3w 2 sin <p t ¡- 2

fr « = Re2j^ , 2 ^ 2 {a1[-(1 + h

1 + h2 + k2

k2) sinL + 2hk cos L] + a2[(1 - hi2 + k2) cos L -

2hksin L] + a3 • 2(hcos L + A-sin L)| •2k -

fj2,n Re2 J2 sin 2 r

a2 • 2h + a3(1 -h2 -k2)]/(1 + h2 + k2)

sin (b =

1 + hi2 + k2

|a1[(1 + h2 -k2)cos Li +

2hksin L ] + a2[(1 - h2 + P)sin L + 2hkcos L] + a3 • 2(hsin L - kxos L)] (28)

5. Multiple-shooting Technique

The multiple-shooting technique is used primarily in order to improve NLP convergence robustness. For the Earth-orbit transfer problems, the multiple-shooting technique falls into two categories: one for thrusting arcs and the other for burn/coast sequences. The former is suitable for minimum-time, single-burn transfers, whereas the latter is suitable for minimum-fuel transfers with multiple burn arcs.

5.1. Multiple shooting for thrusting arcs

In this work, the fact is taken into account that a number of discrete state nodes are inserted during a burn arc integration process. The multiple shooting for thrusting arcs is illustrated in Fig.1, where three state nodes at the time points tj, t2, and t3 for a and e are

equally inserted between t0 and tf . The values of

these nodes must be equal to those obtained by integrating the governing equations at the preceding time point. For example, at each time point t1, t2, and t3,

we need to transform the equinoctial elements computed from numerical integration into classical orbital elements; record the differences between the computed and the guessed values of a and e, then replace a and e by the guessed values while keeping other classical orbit elements the same; transform the classical orbital elements to equinoctial states; continue the integration to the next time point by repeating the process. Finally, the differences between the computed and the guessed values of a and e should be made to zero through the equality constraints in NLP by treating the guessed values as NLP variables. It can be seen that multiple shooting is only applied to two elements, not six. In fact, multiple shooting with part of state nodes might be good enough to solve specific transfer problems. Furthermore, the nodal values might be defined in equinoctial elements so that transforming to and from classical orbit elements (see Ref.[21] for transformation) becomes unnecessary. However, guessing the state nodes as classical orbital elements might be more intuitive for certain types of orbit transfers.

For each burn arc, assuming that n state nodes 55(t1), 55(t2), —, 5c(tn) as initial guesses are inserted between t0 and tf , and the states at time points t1, t2, —, tn, obtained by numerically integrating the equations of motion, are denoted by 5c'(t1), 5'(t2), —,

■5 '(tn).

The extra equality constraints in NLP take the following form:

#ms(t;) = 5; %) - 55 (ti) (i = 1,2,.», n) (29)

For a single burn arc, if q states or classical orbital elements are inserted into the burn duration, there are a total of qxn extra free variables and q*n extra equality constraints in NLP.

State nodes for multiple shooting Control steering nodes ► Integration process

tu t\ ti h tf

Fig.1 Illustration of multiple shooting for thrusting arcs.

5.2. Multiple shooting for burn/coast sequences

Without loss of generosity, the initial guess could be provided for the terminal states of all or a subset of dynamic stages specified in the flight sequence. For Earth-orbit transfers, the main dynamic stages are burn

and coast arcs. The terminal states of each burn or coast arc can be treated as NLP variables, and the guessed values are constrained to the values obtained by computing the corresponding dynamic stage. The guessed terminal states together with the terminal states, which are not optimization variables, are assumed to be the initial condition for the next dynamic stage. Supposing there are m flight sequences, the extra equality constraint is denoted by

Htraj (j) = 5f(j) - 55f( j) (j = 1,2, — , m) (30) where 55f- (j) is the guessed terminal states of the jth dynamic stage (also NLP variables), and 55f( j) the

computed terminal states of the jth dynamic stage. The number of equality constraints depends on the number of dynamic stages and the number of terminal states for each stage. Terminal states can be expressed in classical orbit elements or equinoctial elements.

5.3. Summary of multiple shooting

The main purpose of multiple shooting is to improve convergence robustness without significantly increasing the dimensional size of the NLP problem. By inserting state nodes to constrain the states within proper values, and by properly choosing the terminal states for burn/coast arcs, the entire nominal trajectory does not deviate too much from the optimal solution. In fact, a big terminal error from single shooting is to be replaced by a number of small errors at intermediate time instants from multiple shooting. Moreover, the nodal values might be more intuitive to guess in equinoctial elements than in classical orbital elements. Also, in multiple shooting, a subset of states or all states can be guessed, and the nodes could be equally or unequally spaced.

6. Numerical Examples

The following illustrative examples include transfers from low Earth orbit (LEO) to high elliptic orbit (HEO), from LEO to medium Earth orbit (MEO), from LEO to geostationary Earth orbit (GEO), and from geostationary transfer orbit (GTO) to GEO. Table 1 shows the classical orbit elements for LEO, HEO, MEO, GTO, and GEO. The intermediate thrust acceleration in this article is set to be (T / W )0 on the order of 103, 102, or above. The transfer trajectories usually consist of a few revolutions with relatively large changes in orbit elements. The direct control parameterization and multiple shooting are used to solve the optimal transfer problem and the trajectory optimization is performed in the transformed inertial coordinate. The solutions to be achieved are referred to the J2000 inertial coordinate unless otherwise stated. The minimum-time and minimum-fuel transfers are obtained with the J2 perturbations included unless otherwise stated. The fixed-step fourth-order Runga-Kutta

method is employed to integrate the equations of motion. If the interpolation through the control steering nodes is of linear type, the control steering of a converged solution might be not smooth enough. In this work, the linear interpolation is used to acquire converged solutions at the first step, which are used to be the initial guess in a new optimization problem where the cubic spline interpolation is used to obtain smoother control steering profiles. In all cases the circumferential control defined by Eq.(11) is utilized to be the initial guess for the control steering nodes. The multiple-shooting variables are selected in different ways depending on different transfer problems. In all cases, the initial spacecraft mass is 6 100 kg, and the Isp 3 800 s unless otherwise stated. The final mass obtained in the following examples is also represented by the fraction of the initial mass (final-to-initial mass ratio mf / m0 ).

Table 1 Orbital elements of LEO, HEO, MEO, and GEO

Orbital elements LEO HEO MEO GTO GEO

a/Re 1.047 0 4.076 4 4.076 4 3.820 0 6.610 7

e 0.010 0.700 0 0.731 0

i/(°) 28.5 60.0 Varying, > 60.0 28.5 0

s/(°) 0 30 Free 0 N/A

®/(°) 0 20 N/A Varying N/A

ö/(°) Free Free N/A Free N/A

6.1. Minimum-time, LEO-HEO and LEO-MEO transfers

In this section, the inertial coordinate transformation and multiple shooting for thrusting arcs are demonstrated to effectively solve minimum-time orbit transfer problems. The first example is about transferring the spacecraft from LEO to HEO, which is a typical transfer between two elliptical orbits without allowable coast arcs. With a = Qi (30°), J3 = if (60°), and y = wf (20°), the terminal equinoctial elements are pf = af(1 -ef2) = 2.079 0Re , ff = 0.7 , gf = 0 , and hf = kf = 0 in the transformed inertial coordinate in which the initial classical orbital elements are a0 =

1.047Re, e0 = 0.01:

/q = 37.13°

Í20 = 183.28°

(o0 = 134.16°. In this example, the multiple-shooting

variables are equally spaced along the time axis and generated automatically under the initial and terminal conditions in a linear manner. Only the first five equinoctial elements are used for multiple shooting. The optimal control problem is then converted to an NLP one where a total of 152 variables is optimized, which includes 120 nodes for the control steering, each direction cosine having 40 nodes, 30 nodes for multiple shooting, each of the first five equinoctial elements having six nodes, burn arc duration, and the initial true anomaly of the LEO. Table 2 presents the optimal so-

lutions of LEO-HEO transfers with different initial-to-weight ratios. All the solutions are obtained from the simple initial guesses that are circumferential control steering and reasonable burn arc durations.

In the case with the (T / W)0 of 5x10"3, the time

histories of equinoctial elements at the first NLP iteration in the transformed inertial coordinate are presented in Figs.2-4, in which the nodes for multiple shooting at the first iteration are denoted by circle marks. It can be seen that the multiple-shooting technique drives the equinoctial elements to come close to the target values. In fact, multiple shooting converts one big terminal error into a few small errors that are easier to turn zeros through NLP. Fig.5 presents the optimal trajectory, which consists of about nine revolutions, and Fig.6 shows the control steering.

Table 2 Optimal solutions of LEO-HEO transfers

(T/W)q Transfer time/h m/kg mf/m0 Optimized true anomaly/(°)

10"1 2.01 4 937.90 0.809 5 94.08

5x10"2 3.40 5 117.48 0.838 9 337.88

10"2 15.73 5 190.94 0.850 1 309.76

5x10"3 31.11 5 201.11 0.852 6 321.45

Fig.2 Equinoctial elementp in the transformed inertial coordinate ((T/W)0 = 5x10"3).

0.6 0.4 0.2 0 -0.2

.- With multiple shooting f

----- Without multiple shooting «

r ? g .

Fig.3 Equinoctial elements f and g in the transformed inertial coordinate ((T/W)0 = 5x10"3).

The second example is a series of LEO-MEO transfers with large inclination changes. The (T / W )0 is set

at 102. The MEO orbit with 60° inclination is a representative of global-positioning-system orbit. The transfers to MEO orbits with inclination of 90° and 180° are also investigated. Since the terminal RAAN is free, and the argument of perigee of the target orbit is not

defined, the inertial coordinate transformation requires a = 0, p = if, and y = 0 . The terminal constraints in

the transformed inertial coordinate are pf = af(1 -ef2) = 4.076Re, f = gf = 0 , and hf = kf = 0 . In this

example, the classical orbit elements are employed as the multiple-shooting variables, which require the transformation from classical orbital elements to equinoctial elements (see Ref.[21]). However, only six equally spaced nodes of the elements a, e, and i along the time axis are used for multiple shooting so that the number of NLP variables does not increases too much. As an example, from Table 3, which presents the multiple-shooting variables for transfer to MEO with inclination of 180° in the transformed inertial coordinate, it is noticed that the singularity occurs when the inclination is 180°. In the transformed inertial coordinate, the initial inclination is (180°-28.5°) = 151.5° (2.644 rad), and the final desired inclination is zero. Thus, the singularity does not exist in the transformed inertial coordinate.

Table 3 Multiple-shooting variables for LEO-MEO transfer to the retrograde equatorial orbit

Orbital elements Multiple-shooting variables

a/Re 1.5, 2.5, 3.5, 4.0, 4.5, 4.0

e 0.1, 0.3, 0.5, 0.7, 0.9, 0.1

i/rad 2.6, 2.4, 2.2, 2.0, 1.8, 0.1

Burn time/h 36

Fig.4 Equinoctial elements h and k in the transformed inertial coordinate ((T/W)0 = 5x1Q"3).

: -} I /

«wvy! r..... 7

Fig.6 Control steering for optimal LEO-HEO trajectory ((T/W% = 5x10"3).

Table 4 shows the optimal solutions of LEO-MEO transfers. Figs.7-9 illustrate the time histories of semimajor axis, eccentricity, and inclination. In the cases with 90° and 180° terminal inclinations, the semimajor axis and eccentricity increase so quickly that the inclination near the apogee changes to the desired value more efficiently, and then the semimajor axis and eccentricity are driven back to the desired values. This is a typical feature of transfers with large plane changes. In the case of relatively small plane change (e.g. of 60° terminal inclination), the semimajor axis does not overshoot. Fig.10 shows the optimal trajectory to the retrograde equatorial orbit and Fig. 11 the opti-

Table 4 Optimal solutions of LEO-MEO transfers

Terminal inclination/(°) Transfer time/h mf/kg mf/m0 Optimized true anomaly/(°)

60 14.48 5 262.99 0.862 8 239.17

90 20.36 4 923.36 0.807 1 195.40

180 25.82 4 608.02 0.755 4 100.46

Fig.7 Time histories of semimajor axis for optimal LEO-MEO trajectories.

Fig.5 Optimal LEO-HEO trajectory ((T/W)0 = 5x10 3).

Fig.8 Time histories of eccentricity for optimal LEO-MEO trajectories.

Fig.9 Time histories of inclination for optimal LEO-MEO trajectories.

Fig.10 Optimal LEO-MEO trajectory with final inclination 180°.

Fig.11 Control steering for optimal LEO-MEO trajectory with final inclination 180°.

mal control steering. Although this transfer is not realistic in practices, it stands for the difficulty of a problem that can be solved by the presented direct approach. Actually, the orbit transfers with small plane changes are much easier to solve.

6.2. Minimum-fuel and multiple-burn, LEO-GEO and GTO-GEO transfers

Over decades, many researchers pointed out that the continuous-thrust multiple-burn maneuver might accommodate more payloads than the single-burn maneuver. Different numerical methods were used to attain multiple-burn trajectories, either fully optimal or near-optimal. J. M. Hanson, et al.[13] presented a wide range of multiple-burn trajectories using the TPBVP formulation for each burn arc. J. C. H. Chuang, et al.[14] analyzed the family of multiple-burn transfers using the optimal control theory. Later, they introduced a new method in which the multiple-burn trajectories are

not attained as a whole but, instead, as a sequence of one-burn transfers between optimally chosen intermediate orbits[15]. K. P. Zondervan, et al.[7] utilized a hybrid method to attain 3-burn trajectories with a large plane change. C. A. Kluever, et al.[16] obtained Earth-moon optimal trajectories with up to 5 burns using a direct/indirect approach. They all took advantages of the optimal control theory and TPBVP to solve the problems, that is to say, the costate equations are used to govern the optimal control. In this work, it is demonstrated that the direct approach with multiple shooting is straightforward in solving free or fixed end-time multiple-burn Earth-orbit transfer problems.

First, LEO-GEO trajectories with up to 12 burn arcs are demonstrated. The (T / W )0 is set at 2x10"2 and

the spacecraft mass at GEO is to be maximized. For transfers with a large number of burn arcs, the durations of burn arcs are very short that they might be less than one orbital revolution. Thus, the multiple shooting for burn/coast sequences is more helpful to acquire optimal trajectories than the multiple shooting for thrust arcs. Since GEO is circular and equatorial, the rotating Euler angles are all set to be zeros. As usual, a flight sequence with consecutive burn and coast arcs are specified in advance. At first, the Earth J2 perturbations are not included. The coast arc is computed using the two-body dynamics analytically, and the burn arc is numerically integrated. The SQP design variables include control nodes for each burn arc, time duration for each burn arc, angle displacement for each coast arc, and departure true anomaly at LEO. With multiple shooting for burn/coast sequences, the terminal states of each burn and coast arc are also needed to guess. For the 2- or 3-burn LEO-GEO transfers, it is easier to obtain solutions even without using multiple shooting. In this work, the multiple-shooting technique is used for the transfers with more than 3 burns. Semimajor axis, eccentricity, and inclination for the terminal state of each burn arc should be guessed except the last one. The terminal true anomaly for each coast arc should be guessed as well. However, a problem would be encountered that the burn and coast arcs are prone to be driven to zero by NLP if a large number of burn arcs are assumed, which actually results in fewer-burn solutions. To ease this problem, constraints on the lower and upper bounds of the terminal true anomaly should be imposed for each coast arc. This is beneficial to prevent burn and coast arcs from being driven to be zeros. For example, Table 5 lists the initial guesses for the terminal states of burn and coast arcs for a 12-burn LEO-GEO transfer (with 11 coast arcs), which assumes a structure of a number of perigee burns and a last apogee burn thereby resulting in the optimal solutions from the simply defined initial guesses. Note that the initial control steering is in the circumferential direction, and the durations of burn arcs and coasting angles are simply set to be reasonable values. By using the similar strategy for the transfers with four or more burn arcs, optimal LEO-GEO transfers with up to 12

burn arcs are obtained. Table 6 summarizes the optimal solutions.

Table 5 Initial guesses for terminal states of burn and coast arcs

Terminal classical orbital elements

Segments -

af/Re ef if/rad

1st burn arc 1.5 0.10 0.48

2nd burn arc 2.0 0.15 0.47

3rd burn arc 2.5 0.20 0.46

4th burn arc 3.0 0.25 0.45

5th burn arc 3.5 0.30 0.44

6th burn arc 4.0 0.35 0.43

7th burn arc 4.5 0.40 0.42

8th burn arc 5.0 0.45 0.41

9th burn arc 5.5 0.50 0.40

10th burn arc 6.0 0.55 0.39

11th burn arc 6.0 0.60 0.38

6f(i ) = 6 rad, i = 1,2, — ,10, and %<0f(i) < 2n

Coast arcs 0f(11) = 3 rad, and 0 < i 9f (11) < 2 n

Table 6 LEO-GEO multiple-burn transfers by using direct-multiple-shooting method

J2 not included J included

Number

of burns Transfer time/h mf/kg mf/m0 Transfer time/h mf/kg mf/m0

2 11.67 5 309.45 0.870 4 11.68 5 309.21 0.870 4

3 13.95 5 344.53 0.876 2 13.95 5 344.39 0.876 1

4 13.82 5 361.07 0.878 9 13.81 5 360.87 0.878 8

5 16.35 5 390.95 0.883 8 16.35 5 390.87 0.883 7

6 19.52 5 407.92 0.886 5 19.51 5 407.85 0.886 5

7 22.89 5 414.29 0.887 6 22.84 5 414.19 0.887 6

8 26.82 5 420.70 0.888 6 26.68 5 420.57 0.888 6

9 29.75 5 428.96 0.890 0 29.66 5 428.78 0.890 0

10 33.25 5 432.05 0.890 5 33.28 5 431.85 0.890 5

11 36.78 5 434.30 0.890 9 36.86 5 434.08 0.890 8

12 40.31 5 435.99 0.891 1 40.19 5 435.73 0.891 1

8, 2 apogee burns 36.43 5 416.22 0.887 9 N/A N/A N/A

11, 2 apogee burns 46.65 5 430.14 0.890 2 N/A N/A N/A

The feature of these transfers lies in that the all trajectories consist of certain number of perigee burns and an apogee burn. Table 6 indicates that more burn arcs would result in better performances (for transfers with one apogee burn). However, the time needed by transfers with more burns might be unnecessarily longer than transfers with fewer burns, as pointed out by S. Chuang and R. Goodson[11] (see 3- and 4-burn solutions in Table 6). By loosing the constraints on the terminal true anomaly for each coast arc, 8- and 11-burn transfers with two apogee burns are obtained as shown in Table 6. However, it is found that with the same number of burn arcs, the transfers with two apogee burns would spend longer time and more fuel than those with only one apogee burn for the LEO-GEO transfer problem.

In this case, the burn arcs of the transfers with more than 3 burns are all less than one orbit revolution, whereas the first burn arc of the 2- or 3-burn solution is more than one revolution. The transfer time increases almost linearly once the number of burn arcs exceeds four, but the final mass gets slightly decreasing increments as the number of burn arcs increases. It is noticeable that although solving transfers with more than 12 burns is within the bounds of possibility, very small improvements in performances would be anticipated, also making the assumed burn and coast arcs prone to converge to zeros in NLP problems.

If taking account of J2 perturbations, the coast arcs should be numerically integrated. With solutions that do not consider the J2 perturbations, it is easy to compute the duration of each coast arc. The transfers with the J2 perturbations can be obtained with the solved solutions and the duration of each coast arc as the initial guesses. Table 6 also presents the optimization results. It is understood that the J2 perturbations do not significantly affect the performances of LEO-GEO transfers with the (T / W)0 of 2x10"2. Figs.12-15 present the LEO-GEO transfers with 3, 4, 12, and 8 burn arcs (two apogee burns) respectively. Figs.16-18 show the semimajor axis, eccentricity, and inclination. Figs. 19-20 illustrate the control steering for 3- and 12-burn transfers. It is shown that the perigee burn arcs are almost in line with the velocity direction to efficiently increase the semimajor axis, and the maximum yaw angle is close to the nodal crossing. The eccentricity continuously increases by the perigee burns and reduces to zero only by the apogee burns, and almost all the inclination changes occur during the apogee burns.

Burn Coast

Fig.12 Optimal 3-burn LEO-GEO trajectory (J2 included)

Fig.13 Optimal 4-burn LEO-GEO trajectory(J2 included).

Fig.14 Optimal 12-burn LEO-GEO trajectory (J2 included).

Fig.15 Optimal 8-burn LEO-GEO trajectory with two apogee burn arcs (J2 not included).

Fig.16 Time histories of semimajor axis for LEO-GEO transfers.

Fig.17 Time histories of eccentricity for LEO-GEO transfers.

Fig.18 Time histories of inclination for LEO-GEO transfers.

Fig.19 Control steering for optimal 3-burn LEO-GEO transfers (J2 included).

Fig.20 Control steering for optimal 12-burn LEO-GEO transfers (J2 included).

However, the optimal burn structure would be different if the initial orbit is highly eccentric just as GTO. A GTO-GEO transfer with the (T / W)0 of 5x 10 3 is

taken as an example, and the J2 perturbations are not included in this case. With the multiple shooting for burn/coast sequences and the initial guesses for the terminal semimajor axis, eccentricity, and inclination of each burn arc except the last one and terminal true anomaly of each coast arc, the GTO-GEO transfers with up to 5 burn arcs are solved when a)0 = 0°. Table 7 summarizes the optimal solutions. With up to 5 burn

arcs, the performance improvement is found to be extremely small, and the optimal control structure is a group of apogee-centered burn arcs in this case because the eccentricity vector is assumed to be in line with the nodal crossing line for GTO, thus enabling the large eccentricity and inclination corrections to be most efficiently performed around apogees. However, if the GTO's apogee is not close to the nodal crossing line, a reference to the worst case occurring as (D0 = 90°, 270°, the fuel consumption is not as good as in the case as co0 = 0°, 180° since the eccentricity and inclination can not change most efficiently and simultaneously. Based on the solutions in Table 7, a series of results are acquired with a number of values of between 0° and 90° (see Table 8). It is shown

that the farther the eccentricity vector is away from the nodal crossing line, the worse the performance becomes. It is noteworthy that if the argument of perigee

of GTO are &>0, 180° -&>0, 180° + and 360° - g)0,

the symmetric transfer orbits with the same performance can be obtained. Therefore, it is recommended that the initial elliptic orbit should have the eccentric vector in line with the nodal crossing line if a relatively large plane change is desired. The optimal 5-burn trajectory is presented in Fig.21, and the 5-burn transfer trajectory with = 90° is presented in Fig.22, where the locations of burn arcs are different from those of the transfer trajectory when = 0° . Figs.23-25 show the time histories of semimajor axis, eccentricity, and inclination, respectively.

Table 7 Optimal solutions of GTO-GEO transfers with up to 5 burn arcs

Number of burn Transfer time/h „ , Optimized true mf/kg mf/m0 anomaly/(°)

2 19.21 5 781.05 0.947 7 155.84

3 31.29 5 796.20 0.950 2 165.58

4 44.77 5 800.85 0.951 0 169.60

5 58.84 5 802.90 0.951 3 171.89

Table 8 Optimal solutions of GTO-GEO transfers with different values of w0

Number of ^0=0° ^0=22.5° ^0=45°

burn arcs mf/kg mf/m0 mf/kg mf/m0 mf/kg mf/m0

2 5 781.05 0.947 7 5 743.25 0.941 5 5 690.11 0.932 8

3 5 796.20 0.950 2 5 751.38 0.942 8 5 702.36 0.934 8

4 5 800.85 0.951 0 5 752.69 0.943 1 5 704.68 0.935 2

5 5 802.90 0.951 3 5 753.14 0.943 1 5 712.45 0.936 5

Number of ^0=67.5° ^0=90°

burn arcs mf/kg mf/m0 mf/kg mf/m0

2 5 655.48 0.927 1 5 652.12 0.926 6

3 5 668.55 0.929 3 5 654.54 0.927 0

4 5 670.52 0.929 6 5 655.67 0.927 2

5 5 676.43 0.930 6 5 658.86 0.927 7

Fig.21 Optimal 5-burn GTO-GEO trajectory when ^0= 0° (J2 not included).

Fig.22 Optimal 5-burn GTO-GEO trajectory when <^0= 90° (J2 not included).

Fig.23 Time histories of semimajor axis for GTO-GEO transfers (J2 not included).

Fig.24 Time histories of eccentricity for GTO-GEO transfers (J2 not included).

Fig.25 Time histories of inclination for GTO-GEO transfers (J2 not included).

Multiple burns are also investigated for the fixed-time transfers. It is interesting to determine the optimal number of burn arcs for a fixed-time transfer. After having obtained several solutions with different number of burns for a free-time transfer, the best number of burns for a fixed-time transfer intuitively can be determined. For example, if a 15-hour LEO-GEO transfer is desired, the solution could be expected with either 4 burns or 5 burns in terms of the transfer time in Table 6. With the obtained free-time solutions as the initial guesses, the transfer time can be constrained by NLP and a new optimization problem can be formulated. The optimization results show that the final values of mass at GEO are 5 359.50 kg and 5 386.61 kg for 4- and 5-burn transfers. From this, it is concluded that the best number of burns for the 15-hour transfer is five. This is a simple and efficient way to determine a suitable number of burn arcs for fixed-time transfers without verifying the condition of switching function in the TPBVP.

6.3. Earth-orbit transfers with variable specific impulses

For all transfers in Sections 6.1-6.2, the specific impulse Isp is assumed to be constant. However, for advanced electric propulsion, Isp can be modulated so as to render the thrust magnitude and mass flow rate (see Eq.(7)) changeable thereby resulting in less fuel consumption. Refs. [24-25] pointed out the effects of Isp-modulation on transfers with very-low-continuous-thrust acceleration or large number of orbital revolutions. A simple model of Isp-modulation supposes that the thrust magnitude is inversely proportional to Isp, which itself has lower and upper bounds. For a minimum-time transfer, Isp always works at its lower bound resulting in the largest thrust and most fuel consumption. Therefore, Isp-modulation is generally discussed in connection with certain constraints such as the fixed transfer time, which is greater than the minimum value. Similar to direct control parameterization, the time histories of Isp during each burn arc are interpolated through an appropriate number of nodes, which are

treated as NLP design variables.

For the LEO-GEO transfer presented in Section 6.2, it is assumed that Isp lies between 3 000 s and 3 800 s, which corresponds the initial thrust-to-weight ratio (T / W)0 from 2.53x10"2 to 2.00x10 2. If coasting

arcs are not allowable as in the single-burn transfer, the minimum transfer time is 7.31 h. with Isp at 3 000 s and the final mass of 4 745.05 kg, and the maximum final mass is 5 082.44 kg with Isp at 3 800 s and the transfer time of 8.80 h. A new NLP optimization problem could be formulated, in which the transfer time is set to be from 7.31 h to 8.80 h. and the final mass is to be maximized or the propellant fuel minimized, and the solution with the constant specific impulse can be used as the initial guess. A series of fixed-time solutions are obtained and Fig.26 presents the time histories of Isp. It is shown that Isp always tends to operate at its upper bound for the purpose of saving fuel consumption, and the longer the transfer time is, the longer time Isp operates at its upper bound. Moreover, the fuel consumption decreases as the transfer time increases.

\ / 1 V F' \ 'r \ \ - 1 / \

. \ /1 1 1 Í ' / V ' i i \ i i Í 1 ! • ' ! 1 11- III \ \ \ '

' 1 i ! ' 1 1 i 1 < ' 1 I 1 ■ ¡ '< i i i i i i 1 ! 1 1 M — | 1 — 7.32 h

' » J ' i ! ! ' j i 1 \ — 7.50 h "

' ' i ' ' / í \ —8.00 h "

■ ' i 1 ' '< 1 1 if J / \ — 8.50 h -

f 1 / \ \ 1 \ « 70 h

Fig.26 Time histories of specific impulses.

In need of acquiring a transfer with the final mass greater than 5 082.44 kg or the transfer time longer than 8.80 h, the coasting arcs are required. In fact, coasting arcs implies that Isp operates at an infinite value, which is more efficient than itself operating at a finite upper bound from the angle of saving fuel consumption. For a fixed number of burn arcs without limiting the transfer time, the maximized final mass is achieved with the constant Isp at its upper bound. For 2- and 3-burn LEO-GEO transfers cited in this article (see Fig. 12), the length of the first burn arc is more than a revolution, which implies that modulating Isp during the first burn is helpful for saving fuel consumption (similar to Fig.26) with limited transfer time. However, for transfers with more than 4 burns, each burn is over a small orbit arc within each revolution. The more the number of burn arcs, the shorter the burn arcs within each revolution. These relatively short burn arcs are already located at the best places to carry out efficient changes in orbit elements. Even if Isp is optimized, the optimum value of Isp is still close to the upper bound.

In summary, the multiple-burn scheme with the largest constant Isp affords a better choice than /^-modulation if only the purpose of fuel-saving is pursued. The increased transfer time resulting from larger number of burn arcs is on the order of hours for Earth-orbit transfers in this article. In fact, Isp-modulation is more suitable for interplanetary trajectories, which usually last months and years, as such limiting the transfer time becomes much more necessary. In addition, interplanetary trajectories may include many intermediate constraints such as gravity assists, intermediate and rendezvous, which exclude totally free transfer time.

7. Comparisons with Other Methods

Compared with the collocation method proposed by A. L. Herman and D. B. Spencer[5], the direct control parameterization usually requires fewer NLP variables and constraints since the states are not parameterized. Another difference lies in that the direct control parameterization uses the explicit integration method, whereas the collocation method uses certain curve fitting methods to approximate the time histories of states. The solution accuracy obtained with the collocation method usually needs verification afterwards.

Y. Gao, et al.[8,26] presented a hybrid method in which both state and costate equations are integrated with multiple shooting. In this method, the initial co-state variables for each burn arc need to be guessed. In general, the convergence robustness of the hybrid method lies between the direct and indirect methods. For a single-burn transfer, the NLP design variables can be reduced by using the hybrid method in place of the direct method, and this advantage turns more significant as the transfer time becomes longer[26]. For multiple-burn transfers in this study, each burn arc is relatively short and the control steering appears almost linear. Therefore, two or three nodes are needed to interpolate into the continuous control steering. In this case, the hybrid method does not possess the advantage to reduce number of design variables compared with the direct control parameterization.

Compared with nongradient-based optimization methods or stochastic methods, such as genetic algorithm^7-1 and differential evolution[18], the most outstanding feature is that the gradient-based direct method only finds local optima, whereas the stochastic method finds global ones. However, on the flip side, the stochastic method usually requires longer searching time and lacks absolute certainty of finding out global optima. Consequently, it is usually suggested to use the stochastic method for coarse searching and the gradient-based optimization for finding accurate solutions. For orbit transfers consisting of only burn and coast arcs under study, the gradient-based optimization is more efficient than the stochastic methods. As for the very complex interplanetary transfers inclusive of a variety of system parameters, the stochastic methods

are beneficial for finding coarse solutions which are optimized further with the gradient-based methods.

8. Potential Guidance Scheme Using Model Predictive Control

Model predictive control (MPC), also called receding horizon control, is meant to repeat computation of the control inputs by solving the optimal control problem over a given horizon, thus producing an open-loop optimal control sequence. The first control over the total time span in the sequence is applied at the outset. In the next sampling instant, a new optimal control problem arises and is to be solved based on new measurements. MPC has found successful applications in many areas such as the process control industry. The main drawback of MPC is the prolonged repeated calculation of new optimal controls for the real-time control system. However, this drawback turns less significant for orbit transfers, which always spend hours or even days. Another key problem of MPC is about the stability or the optimization convergence of NLP for trajectory optimization. The proposed direct method has illustrated good convergence robustness for a wide range of transfer problems with simply defined initial guesses, which provides MPC with a promising prospect in the autonomous guidance of spacecraft.

The preliminary implementation of the autonomous guidance scheme by way of MPC could be summarized as follows:

(1) To compute an optimal control offline to obtain optimal solutions of control nodes and state nodes for multiple shooting. The optimal parameters are up-linked to the spacecraft and used for real control during the first time horizon.

(2) In the next sampling instant, to re-allocate control nodes and state nodes on the base of prior-obtained solutions and to re-optimize control steering onboard by using the re-allocated control nodes and state nodes as the initial guesses.

(3) To apply the new solution for the optimal control during the succeeding time horizon, then return back to Step (2) and repeat the same process until the target orbit is reached.

Except for the offline optimization at the first step, all the optimizations in the Step (2) are performed with good initial guesses. Thus, NLP optimization should converge quickly and robustly. Of course, the number of control nodes and state nodes, which affects the NLP computation time, should be appropriately chosen for the onboard optimization. The detail implementation and verification of the autonomous guidance scheme using MPC is reserved for future research.

9. Conclusions

The proposed direct approach is particularly suitable to solve a variety of Earth-orbit minimum-time and

minimum-fuel transfers with intermediate acceleration. The inertial coordinate transformation strategy is employed to produce general simple terminal constraints and remove the singularity occurring at the retrograde equatorial orbit. Moreover, the simple terminal constraints also improve the convergence robustness of nonlinear programming. The strategy to deploy the multiple-shooting variables is described for Earth-orbit transfer problems. Numerical examples in the article include minimum-time, minimum-fuel trajectories with up to 12 burn arcs, as well as fixed-time transfers with variable specific impulses. All the solutions are obtained with multiple shooting in the transformed iner-tial coordinate from simply defined initial guesses. As for fuel-saving Earth-orbit transfers, it is concluded that the multiple-burn scheme with largest constant Isp offers a better choice than Isp-modulation.

The proposed direct approach is also easy to solve the transfers with lower acceleration. However, trajectories with hundreds or thousands of revolutions require a large number of nodes for control steering and steps for integrating the equations of motion, which results in heavy computational loads. For transfers with intermediate-level acceleration under study, shorter computation time and improved convergence make it possible to use the model predictive control in the onboard guidance scheme.

References

[1] Rayman M D, Williams S N. Design of the first interplanetary solar electric propulsion mission. Journal of Spacecraft and Rockets 2002; 39(4): 589-595.

[2] Rathsman P, Kugelberg J, Bodin P, et al. SMART-1: development and lessons learnt. Acta Astronautica 2005; 57(2-8): 455-468.

[3] Kechichian J A. Optimal low thrust orbit geostationary Earth-orbit intermediate acceleration orbit transfer. Journal of Guidance, Control, and Dynamics 1997; 20(4): 803-811.

[4] Herman A L, Conway B A. Direct optimization using collocation based on high-order-gauss-lobatto quadratic rules. Journal of Guidance, Control, and Dynamics 1994; 17(3): 480-487.

[5] Herman A L, Spencer D B. Optimal low-thrust earth-orbit transfers using higher-order collocation methods. Journal of Guidance, Control, and Dynamics 2002; 25(1): 40-47.

[6] Scheel W A, Conway B A. Optimization of very-low-thrust, many-revolution spacecraft trajectories. Journal of Guidance, Control, and Dynamics 1994; 17(6): 1185-1192.

[7] Zondervan K P, Wood L J, Caughey T K. Optimal low-thrust, three-burn orbit transfers with large plane changes. Journal of the Astronautical Sciences 1984; 32(4): 407-427.

[8] Gao Y, Kluever C A. Low-thrust interplanetary orbit transfers using hybrid trajectory optimization method with multiple shooting. AIAA-2004-5088, 2004.

[9] Ilgen M R. Hybrid method for computing optimal low thrust OTV trajectories. Advances in the Astronautical Sciences 1994; 87(2): 941-958.

[10] Kluever C A, Oleson S R. Direct approach for computing near-optimal low-thrust earth-orbit transfers. Journal of Spacecraft and Rockets 1998; 35(4): 509515.

[11] Geffroy S, Epenoy R. Optimal low-thrust transfers with constraints—generalization of averaging. Acta Astronautica 1997; 41(3): 133-149.

[12] Gao Y. Near-optimal very low-thrust earth-orbit transfers and guidance schemes. Journal of Guidance, Control, and Dynamics 2007; 30(2): 529-539.

[13] Hanson J M, Duckman G A. Optimization of many-burn orbital transfers. Journal of the Astronautical Sciences 1997; 45(1): 1-40.

[14] Chuang J C H, Goodson T D, Hanson J. Multiple-burn families of optimal low-and-medium-thrust orbit transfer. Journal of Spacecraft and Rocket 1999; 36(6): 866874.

[15] Goodson T, Chuang, J C H, Hanson J. Optimal finite thrust orbit transfers with large number of burns. Journal of Guidance, Control, and Dynamics 1999; 22(1): 139-148.

[16] Kluever C A, Pierson B L. Optimal low-thrust earth-moon transfers with a switching function structure. Journal of the Astronautical Sciences 1994; 42(3): 269283.

[17] Kim Y H. Spencer D B. Optimal spacecraft rendezvous using genetic algorithms. Journal of Spacecraft and Rockets 2002; 39(6): 859-865.

[18] Dachwald B. Optimization of very-low-thrust trajectories using evolutionary neurocontrol. Acta Astronautica 2005; 57 (2-8): 175-185.

[19] Hull D G. Conversion of optimal control problems into parameter optimization problems. Journal of Guidance, Control, and Dynamics 1997; 20(1): 57-60.

[20] Walker M J H, Ireland B, Owens J. A set of modified equinoctial orbit elements. Celestial Mechanics 1985; 36(4): 409-419.

[21] Betts J T. Optimal interplanetary orbit transfers by direct transcription. Journal of the Astronautical Sciences 1994; 42(3): 247-268.

[22] Pouliot M R. CONOPT2: a rapidly convergent constrained trajectory optimization program for TRAJEX. GDC-SP-82-008, 1982.

[23] Kluever C A. Optimal low-thrust interplanetary trajectories by direct method techniques. Journal of the As-tronautical Sciences 1997; 45(3): 247-262.

[24] Kluever C A. Geostationary orbit transfers using solar electric propulsion with specific impulse modulation. Journal of Spacecraft and Rockets 2004; 41(3): 461466.

[25] Seywald H, Roithmayr C M, Troutman P A, el al. Fuel-optimal transfers between coplanar circular orbits using variable-specific-impulse engines. Journal of Guidance, Control, and Dynamics 2005; 28(4): 795800.

[26] Gao Y. Advances in low-thrust trajectory optimization and flight mechanics. PhD thesis. University of Missouri-Columbia, 2003.