Available online at www.sciencedirect.com Chinese

ScienceDirect leronautfcs

Chinese Journal of Aeronautics 22(2009) 426-433 www elsevier com/locate/cja

Direct Optimization of Low-thrust Many-revolution Earth-orbit Transfers

Gao Yang*

Academy of Opto-Electronics, Chinese Academy of Sciences, Beijing 100190, China Received 22 July 2008; accepted 10 February 2009

Abstract

Low-thrust Earth-orbit transfers with 10~5-order thrust-to-weight ratios involve a large number of orbital revolutions which poses a real challenge to trajectory optimization. This article develops a direct method to optimize minimum-time low-thrust many-revolution Earth-orbit transfers. A parameterized control law in each orbit, in the form of the true optimal control, is proposed, and the time history of the parameters governing the control law is interpolated through a finite number of nodal values. The orbital averaging method is used to significantly reduce the computational workload and the trajectory optimization is conducted based on the orbital averaging dynamics expressed by nonsingular equinoctial elements. Furthermore, Earth's shadowing and perturbation effects are taken into account. The optimal transfer problem is thus converted to the parameter optimization problem that can be solved by nonlinear programming. Taking advantage of the mapping between the parameterized control law and the Lyapunov control law, a technique is proposed to acquire good initial guesses for optimization variables, which results in enlarged convergence domain of the direct optimization method. Numerical examples of optimal Earth-orbit transfers are presented.

Keywords: low thrust; orbital transfer; optimization; direct method; orbital averaging; Lyapunov method

1. Introduction

Compared to vehicles propelled by conventional chemical propellant, spacecraft driven by low-thrust electric propulsion can deliver more payload fraction due to its high specific impulse. However, low-thrust propulsion usually results in long-duration Earth-orbit transfers involving hundreds or even thousands of orbital revolutions. Because of this, optimization of low-thrust transfer trajectories has long been remaining a real challenge. Furthermore, the optimization problem would be further complicated by shadowing and perturbation effects that should be taken into consideration if the spacecraft transfers through low-altitude orbital space.

In the early 1960's, T. N. Edelbaum[1] presented a well-known analytic solution to estimate flight time and total impulse for low-thrust circle-to-circle orbital transfers. SEPSPOT[2], a program developed for Earth-orbit transfers by using solar electric propulsion,

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

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

Foundation item: National Natural Science Foundation of China (10603005)

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

solved the two-point boundary-value problem (2PBVP) with a shooting method. However, the solution to the 2PBVP is sensitive to the initial guess for the costate variables that do not imply any intuitive physical meanings. Compared to the approach of solving the 2PBVP (traditionally termed indirect method), direct methods generally exhibit a larger radius of convergence domain. J. T. Betts[3] used the direct collocation method to solve a 578-revolution transfer problem, thereby bringing forth an optimization problem with 416 123 variables and 249 674 constraints. J. T. Betts' work leads to large-scale nonlinear optimization problems that require tremendous computational workload. W. A. Scheel, et al.[4] developed a parallel Runge-Kutta method for solving many-revolution orbital transfer problems. However, the parameterization of the control steering program in their method still requires a large number of nodes which are used as optimization variables.

In order to alleviate computational workload, both indirect and direct methods adopt orbital averaging used in SEPSPOT. C. A. Kluever, et al.[5] put forward a method that combines orbital averaging and a blended control strategy to solve minimum-time Earth-orbit transfers. S. Geffroy, et al.[6] solved orbital transfer problems with constraints by using orbital averaging.

Recently, Y. Gao[7] developed a new control strategy which also uses the orbital averaging method to solve near-optimal minimum-time and fuel-efficient orbital transfers in the presence of J2 perturbations and shadowing.

In addition to traditional indirect and direct methods, there exist other methods to solve low-thrust many-revolution orbital transfers. For instance, the Lyapunov control can be used to guide the spacecraft to the target orbits[8-10]. However, the shadowing was neglected in these studies. Besides, the Lyapunov control law is usually not optimal. T. Haberkorn, et al.[11] developed a homotopic approach to solve fixed-time minimum-fuel orbital transfers in an ideal gravity field without considering shadowing. The homotopic approach might take a few hours execution time to obtain optimal solutions.

This article develops a direct optimization method to solve minimum-time low-thrust many-revolution Earth-orbit transfer problems. A parameterized control law in the form of true optimal control is formulated in each orbit of the transfer trajectory. The orbital averaging method is used to ease the computation for trajectory propagation. Moreover, perturbations and shadowing can be easily modeled by orbital averaging. The parameters governing the control steering program in each orbit are interpolated through a finite number of nodal values in terms of the time history. The optimal control problem is then converted to the parameter optimization problem that can be solved by nonlinear programming (NLP). In addition, multiple shooting of mean orbital elements is employed to improve the convergence robustness of NLP. Furthermore, the initial guess for the optimization variables of NLP is obtained using a mapping between the parameterized control law and the Lyapunov control law.

2. System Dynamics Based on Orbital Averaging

A set of nonsingular equinoctial elements x = [p f g h k L]T is utilized to govern the equations of motion[12] where p is the semi-latus rectum, L is the true longitude, and the other four do not have intuitive physical meanings. The equinoctial elements can be obtained in terms of the classical orbital elements with p = a (1-e2), f= ecos(®+^), g = esin(®+^), h = tan(i/2)cos Q k = tan(i/2)sin Q and L = f2+a+d where a is the semi-major axis, e the eccentricity, i the inclination, Q the longitude of ascending node, œ the argument of perigee, and 0the true anomaly.

The Gauss variational equation in a matrix form can be written as

x=m ^ ma+fp j+d (1)

where M is a 6x3 matrix and D a 6x1 vector, T the thrust magnitude, m the spacecraft mass, a the unit vector in thrust direction in the local-vertical-local-horizon frame, and fp the perturbation-vector. If the

perturbation is conservative, another equation of motion, i.e. the Lagrange planetary equation[12] holds. In addition, the mass flow rate is

m = -T /( ge IJ

where ge (ge = 9.806 65 m/s ) is the gravitational acceleration at the sea level and Isp the specific impulse. Note that the thrust magnitude T = (T/W)0m0ge is determined by the initial thrust-to-weight ratio (T/W)0 and the initial spacecraft mass m0.

Assuming that the low-thrust acceleration is much lower than the gravitational acceleration, the following can be obtained

-T-ylMP

1 + f cos L + g sin L

where ¿u (¿u = 398 601 km /s ) is the Earth gravitational parameter. The first-order time rates of averaged incremental changes for the first five equinoctial elements (defined as x = [p f g h k ]T, termed mean equinoctial elements) can be computed by

dx 1 f2^ dt ,T

-I x—dL (4)

x—dL dL

where P is the orbital period. The time rates of the mean equinoctial elements are obtained by expanding Eq.(4):

— = x =— (1 - f2 -dt l 2k j '

C l2 f/( x, L,a)

(1 + f cos L + g sin L)

r2)3/2 •

fiP( X, L, fv)

(1 + f cos L + g sin L)2

where x = [x1 x2 x3 x4 x5]T = [p f g h k]T and l =1,2, — ,5. Likewise, the averaged mass flow rate is

dm ^ 1 72 . _ = = __(1 -f -,

Jil (1 + f cos L + g sin L)2

-2\3/2 ' ) '

where fj( x, L,a)and flP( x, L, /p)are the dynamics of the osculating elements due to thrust and perturbation, respectively. The termfj (x,L,a) is computed through the Gauss variational equation (T / m)Ma, where M is a 5^3 matrix after removing the 6th row of M. The term flP( x, L, /p)can be computed through either the Gauss variational equation or the Lagrange planetary equation.

The integral limits L1 and L2 in Eqs.(5)-(6) represent the starting and ending angles (measured in true longitude) of a burn arc in each orbital revolution. By taking the Earth shadow into consideration (null thrust in shadow), the exit and entry angles are logically denoted by L1 and L2, respectively (see Fig.1). Note that, compared with completely continuous thrust, the dis-

continuous thrust caused by the shadow may strongly affect the time evolution of orbital elements. The definite finite integrals in Eqs.(5)-(6) can be approximately computed by the Gauss-Legendre quadrature. Eqs.(5)-(6) can be integrated with large time interval (e.g. one or two days). Compared with propagating osculating elements, propagating mean equinoctial elements considerably lessen the computational workload. For the details of orbital averaging, refer to Refs.[13]-[14].

Vernal equinox dlection Fig.1 Illustration of Earth shadow exit and entry angles.

3. Parameterized Control Law and Nonlinear Parameter Optimization

According to the Calculus of variations or the optimal control theory[15], with a cost function in the Mayer form, the Hamiltonian for the averaged dynamic system (Eqs.(5)-(6)) is defined as

H = XT x + Xmm

where X =[ Ap Af Ag Xh Xk ]T and Am are the

costate variables associated with the corresponding mean equinoctial elements and the spacecraft mass, respectively. The optimum value of the thrust direction a is obtained by setting dH / da = 0 with the constraint aTa=1.

= _ m Тл

a°pt " Il M TX||

In this study, it is assumed that the optimal control steering is directly governed by the parameters X = [ lp If Ag lh \ ]T, and the time history of X(t) is

interpolated through an appropriate number of nodal values along the time axis rather than governed by the costate differential equation (see Ref.[2]). The nodal values for interpolating X(t) should be chosen to be optimal. Therefore, Eq.(8) denotes the parameterized control steering used in each transfer oribt, which is termed the calculus of variations (COV)-based control law. The parameterized COV-based control law possesses a concise and optimal form derived from COV. Note that, different from the case in the traditional

indirect method, the costate differential equations are not used to govern X (t).

The orbital transfer problem is now converted to the parameter optimization problem. With N nodal values (X(tj) (I = 1,2, — , N)) for interpolating the time history of X(t) and tf to denote the transfer time, the parameter optimization problem could be described as follows:

Find optimal variables

[I(ti) I(t2) ••• I(tN) tf ]

to minimize the terminal time tf, subjecting to

® Eqs.(5)-(6);

® x(tf) equals the orbital elements of the target orbit.

It is worth mentioning that the constraint of the equations of motion is automatically satisfied by the explicit numerical integration of Eqs.(5)-(6) from the initial time to the terminal time. In addition, any other optimization variables as well as equality and inequality constraints could be included in above-stated optimization problem. The proposed parameter optimization problem could be solved by NLP—sequential quadratic programming (SQP)[16] .

In order to further improve the convergence robustness of NLP, the multiple-shooting technique has found wide applications in trajectory optimization problems. This article introduces multiple shooting of mean orbital element x . The state nodes xms(tj)

(J=1,2,—,M) at tJ can be equally or unequally inserted between the initial time t0 and the terminal time tf. The values of these nodes must be equal to the values xint(tJ) obtained by integrating the governing equations from the preceding instants. With the multiple-shooting technique, the extra optimization variables are Xms (ti Ь Xms (t2) , Xms (tM ) , and the additional equality constraint is

xms (tj) - Xmt(tj ) = 0 (9)

The state nodes xms(tJ) could be expressed by either classical orbital elements or equinoctial elements. However, guessing the state nodes expressed by classical orbital elements might be more intuitive for some types of orbital transfers. The details for converting an optimal control problem into a parameter optimization problem as well as for using the multiple-shooting technique can be found in Ref.[17].

4. Acquisition of Initial Guesses by Using the Lyapunov Control

Despite the proposed approach being categorized as a direct method, proper initial guesses for optimization variables are always helpful and needed. This section will associate the COV-based control law with the Lyapunov control law firstly made public in Ref.[8] so as to provide a tool to obtain proper initial guesses for

optimization variables.

According to Ref.[8], the Lyapunov function defined by the equinoctial elements can be written into

V = ^[Qp (P - P*)2 + Qf (f - f *)2 +

Qg (g-g*)2 + Qh(h -h*)2 + Qk(k -k*)2]

where the equinoctial elements of the target orbit are denoted by x* = [p * f * g * h * k*]T , and Qp, Q f, Qg, Qh, and Qk are constant values called Lyapunov gains. Now the control law or the Lyapunov control law is defined as

II MTV— |

where V- = Q(x - x*) and Q is a diagonal 5x5 matrix with the diagonal elements Qp, Qf, Qg, Qh, and Qk. Then V could be readily acquired without the consideration of perturbation.

V = VT x=2-(1 - f2 - g 2)3/2 • 2%

rL -(T/m)\\M% |

Jl' (1 + f cos L + g sin L)2

If the gain matrix Q is positive definite, Eq.(12) would imply that the control system involving the control law of Eq.(11) is asymptotically stable (for details see Ref.[8]). To simplify the Lyapunov control law, Q could be defined as a constant time-independent matrix. Because of feedback mechanism, the Lyapunov control law naturally drives the spacecraft to the target orbit. However, the solution obtained through the Lyapunov control law is generally not optimal.

It can be disclosed that the Lyapunov control law and the COV-based control law of Eq.(8) become identical if the following mapping holds true:

%x = Q( x - x*) (13)

Let the guess of Lyapunov gains be QLguess and the transfer time be tguess, a nominal trajectory is propagated with the Lyapunov control law to obtain xL(t). It should be ensured that this trajectory reaches the vicinity of the target orbit, that is to say, xL (tf) gets close to x *. This might require making a few trial-and-errors for determining QL"ess and tguess. With the solutions of Qfessand xL(t), the nodes Xguess(t7) for interpolating the time history of X(t) can be obtained by using the mapping:

Iguess(t7) = QLuess (xL (tj) - x*) (I = 1,2,-, N) (14)

If the multiple-shooting technique is employed, the state nodes for multiply shooting x^ess(tj) can be obtained through xL (t):

Therefore, the initial guesses for the optimization variables in the parameter optimization problem are Iguess(tI),xmT(tJ), and_tfguess. This proposed technique converts guessing Aguess(tI) and xmguess(tJ) into guessing QLuess, which decreases the number of optimization variables to be determined to only five and makes it easier to get the candidate solutions close to the target orbit due to the feedback mechanism of the Lyapunov control.

5. Numerical Examples

This section will present Earth-orbit transfer examples with four types of orbits involved: low Earth orbit (LEO), geostationary transfer orbit (GTO), geostationary orbit (GEO), and high elliptical orbit (HEO). Table 1 lists the classical orbital elements of each orbit where Re is Earth radius and Table 2 three cases of orbital transfers including LEO-GEO, GTO-GEO, and LEO-HEO. Note that both LEO-GEO and GTO-GEO transfers are cited from Refs.[5],[7]. Table 2 also shows the specific impulses, initial spacecraft mass, and initial thrust-to-weight ratios. All the proposed orbital transfers contain the zonal harmonics J2-J5 perturbations and cylindrical Earth shadow. The initial departure date is set to January 1, 2008. The Jet Propulsion Laboratory (JPL) planetary ephemeris[18] is used to obtain the sun's position in computing Earth shadowing. The standard fourth-order fixed step-size Runge-Kutta method is employed to propagate mean equinoctial elements. Two hundred integration steps are set up for propagating the LEO-GEO transfer trajectories and 100 steps for propagating the LEO-HEO and GTO-GEO transfer trajectories.

Table 1 Classical orbital elements of LEO, HEO, GTO, and GEO

Orbit a/Re e i/o m°) o^c)

LEO 1.086 0 0 28.5 0 0

GTO 3.820 0 0.731 27.0 99 0

GEO 6.610 7 0 0 — —

HEO 4.076 0 0.700 60.0 30 20

Table 2 Transfer cases and spacecraft parameters

Transfer case Isp/s mc/kg (T/W)0

LEO-GEO 3 300 1 200 3.413 5x10~5

GTO-GEO 3 300 450 4.551 4x10~5

LEO-HEO 3 300 1 000 8.000 0x10~5

s(tj) = —l (tj) (J = 1,2,-, M)

Just as Section 4 has indicated, the initial guess for the optimization variables can be obtained with the Lyapunov control, which will be demonstrated by taking the LEO-GEO transfer as an example. After a few trial-and-errors, a set of Lyapunov gains—[Qp Qf Qg Qh Qk] = [1 5 5 102 102] is chosen and a number of values for the transfer time are selected. Table 3 lists the terminal errors in semi-major axis, eccentricity, and

inclination (denoted by af - a *, ef - e *, if - i *, respectively) for the four transfer time solutions.

Table 3 Terminal orbital errors with different transfer time

Transfer time/day (af - a *)/Re ef - e * (if - i *)/(0)

200 7.445 0x10"3 3.315 8x10"2 8.437 4

210 8.430 1x10"4 1.159 5x10"2 4.382 2

220 1.443 6x10"2 2.028 5x10"3 0.861 2

230 1.329 2x10"2 1.234 8x10"3 0.153 2

With the simple constant Lyapunov gains, it is observed that the 200-day solution is inferior to the other three ones, which could be used to achieve the initial guesses for X (t) and the multiple-shooting variables.

It should be pointed out that applying the Lyapunov control is not meant to generate an accurate transfer trajectory, but an approximate solution to aquire the initial guesses through the mapping of Eqs.(14)-(15) for further nonlinear optimization problem. In this aspect, with the help of feedback mechanism, rough selection of the Lyapunov gains and the transfer time is capable of achieving favorable approximate solutions.

In the 210-day transfer solution, ten nodes are equally spaced along the time axis for interpolating X(t). The multiple-shooting technique is then used to divide some state variables (semi-major axis, eccentricity, and inclination) into five segments and four nodes are inserted between the initial and terminal time. It is necessary to transform equinoctial elements to classical orbital elements at each multiple-shooting time instant. Therefore, there are a total of 63 SQP design variables, including 50 nodes for interpolating mean costate variables (10 nodes for each element of X(t)), 12 nodes for multiple shooting (four nodes for

each of semi-major axis, eccentricity, and inclination), as well as the transfer time. The equality constraints are five equinoctial elements of the GEO orbit (p * =

6.610 7Re and f * = g* = h* = k * = 0) and 12 constraints resulting from multiple shooting. The same is true of the GTO-GEO and LEO-HEO transfers except for the numbers of optimization variables and constraints which might be different.

In the GTO-GEO transfer, the converged solution could be attained without multiple shooting. The initial guess for optimization variables are generated by the Lyapunov gains [1 10 10 102 102] and the transfer time 75 days. Ten nodes are used for interpolating the time history of X (t) in the optimization problem.

By solving the modeled parameter optimization problems, the converged solutions are obtained. Table 4 presents the optimal solutions of the three transfers. This table also contains the solutions in Ref.[7] (values after "/") for comparison and the approximate orbital

revolutions of each transfer. The time for the LEO-GEO transfer is 199.5 days and for GTO-GEO 66.8 days. In comparison with other previously obtained solutions, they are 198.9 days and 66.6 days by SEPSPOT[5], 200.3 days and 67.1 days by the blended control law[5] and 202.9 days and 70.2 days by the near-optimal control strategy[7]. Figs.2-4 illustrate the time histories of semi-major axis, eccentricity, and inclination for the LEO-GEO transfer and Figs.5-7 for the GTO-GEO transfer. The solutions obtained from the Lyapunov control and the near-optimal solutions in Ref.[7] are also plotted in these figures. In Figs.2-4, the initial values of the state nodes for multiple shooting are marked by circles.

Because the control law used in Ref.[7] is near optimal and the COV-based control law is in the form of true optimal control, the solutions attained in this study are superior to those in Ref.[7] for minimum-time transfers in terms of transfer time and propellant mass. Ref.[7] has evidenced that the very low-thrust transfers

Table 4 Optimal solutions of LEO-GEO, GTO-GEO, and LEO-HEO transfers

Transfer case Transfer time/day Propellant mass/kg Approximate revolutions

LEO-GEO 199.5 / 202.9[7] 193.5 / 197.5[7] 1 249

GTO-GEO 66.8 / 70.2[7] 34.77 / 36.5[7] 96

LEO-HEO 101.0 201.14 913

Fig.2 Time histories of semi-major axis for the LEO-GEO transfers.

Fig.3 Time histories of eccentricity for the LEO-GEO transfers.

Fig.4 Time histories of inclination for the LEO-GEO transfers.

Fig.5 Time histories of semi-major axis for the GTO-GEO transfers.

Fig.6 Time histories of eccentricity for the GTO-GEO transfers.

play minor part in cutting down on propellant consumption (under 19%) but spend almost twice the transfer time reaching as long as 400 days. Therefore, to strike the balance between the transfer time and propellant mass, a minimum-time solution for the very low-thrust orbital transfer might fit the bill.

For the LEO-HEO transfer, the Lyapunov gains are selected to be [1 1 1 20 20] and the transfer time is set to be 110 day. In the optimization problem, the multiple-shooting technique is used only for semi-major axis and nine nodes are equally spaced between the initial and terminal time. Note that as the COV-based control law does not possess feedback mechanism, when the transfer time is much longer than the optimal, the semi-major axis becomes liable to overshoot without using multiple shooting. Figs.8-9 demonstrate the time histories of equinoctial elements for the optimal solution of the LEO- HEO transfer. Note that the equinoctial elements of the target orbit are [2.08tfe 0.45 0.54 0.5 0.29].

Time histories of equinoctial element p for the LEO-HEO transfers.

1 ■ 1 1 1 JK

■ V>' i (V / —/ ■ -g ■ — h -»» k

60 tiday

100 120

Fig.7 Time histories of inclination for the GTO-GEO transfers.

Fig.9 Time histories of equinoctial elementsf g, h, and k for the LEO-HEO transfers.

The thrust direction in each transfer orbit is governed by the parameterized control law defined by Eq.(8). During each integration time interval, i.e. one or two days, the control law in each orbit remains the same. Fig.10 shows the thrust direction angles for the GTO-GEO transfer with the pitch S and yaw ^ defined by

[a1 a2 a3]T = [sin Scos <f> cos Scos <f> sin (16)

where S is measured from the local horizon to the projection of the thrust vector onto the orbit plane and ^ from the orbit plane to the thrust vector. It is indicated in Fig.10 that, from the 30th day to the terminal time of the transfer, the control steering is mainly responsible for circularizing the orbit. In addition, the osculating orbital elements are propagated with the optimization solution to achieve the GTO-GEO transfer trajectory, which is an ascending spiral involving almost a hundred revolutions (see Fig.11, where X, Y, and Z are the three coordinate axes).

Fig.10 Thrust direction angles for the GTO-GEO transfer.

turbed converged solution is used as an initial guess for a new NLP optimization run, a slightly better solution may be obtained. Moreover, a slightly perturbed solution failing to satisfy the convergence criterion is used as an initial guess to restart a new run, which is also beneficial for obtaining converged solutions. This suggests that the designer be required to run optimization a few more times by perturbing existing feasible solutions. Table 5 lists the iteration steps, minimized objectives, and equality constraint norms for the three numerical examples. In the first optimization, the initial guess is provided by trajectory solutions from the Lyapunov control. It takes less than 5 min for implementing the optimization on a PC with a built-in 2.66 GHz CPU.

Table 5 Iterations for the first optimization

Transfer case Iteration step Minimized objective tf /day Equality constraint norm (non-dimensional unit)

LEO-GEO 43 201.05 Converged, 2.794 8x10"7

GTO-GEO 53 66.83 Converged, 1.755 2x10"7

LEO-HEO 57 103.75 Not converged, 2.66 x10"4

For the LEO-GEO and LEO-HEO transfers under study, the solutions obtained in the first optimization do not guarantee the 199.5-day and 101.0-day transfer solutions listed in Table 4. With the above-mentioned suggestions, better solutions could be found after implementing optimization for several times. On obtaining the first converged solution, multiple shooting is usually not needed for subsequent optimization thereby reducing the number of optimization variables. It should be pointed out that different NLP solvers might have different convergence properties.

Fig. 11 GTO-GEO transfer traj ectory.

7. Comparisons with Other Methods

6. Some Comments on Optimization Convergence

Although the direct method generally exhibits a larger convergence domain than the indirect method, simply selected values for optimization variables usually do not always give rise to converged solutions. Note that the LEO-GEO transfer involves much more orbital revolutions (over 1 000) than the other two transfers. With the aid from the Lyapunov control to provide the initial guesses, the practical experience to settle the optimization problem with SQP[16] has shown the possibility of obtaining converged solutions (maybe local optima) in almost all the cases. Otherwise, the optimization of the LEO-GEO is very difficult to converge even with plenty of trials of arbitrarily selected initial guesses for optimization variables X(tj) (I = 1,2,-, N).

All converged solutions obtained by using the direct method are deemed locally optimal. If a slightly per-

In general, the direct optimization method has a larger convergence domain than the indirect method, especially benefited from the initial guess generated by using the Lyapunov control. This advantage is witnessed by comparing with solving costate differential equations in Refs.[2],[6]. With orbital averaging, the computational workload can be enormously reduced if compared with the work in Refs.[3],[4], where direct control parameterization inextricably leads to large dimensional nonlinear optimization problems. Furthermore, the perturbations and shadowing can be easily taken into consideration by using orbital averaging. The parameterized control law is truly an optimal control law as distinct from the Lyapunov control law[8-10], blended control[5], and near-optimal control law[7] that could not be regarded optimal theoretically. Finally, it is worth noting that the use of equinoctial elements eliminates the singularity when either eccentricity or inclination equals zero, which is not taken into account

by Refs.[5],[7] where the control laws with intuitive physical meaning are all expressed in classical orbital elements. The proposed method is effective and fit to deal with the minimum-time transfers. It is anticipated to be extended to tackle the fuel-efficient transfers with active coasting mechanism.

8. Conclusions

A direct approach has been developed to optimize minimum-time low-thrust many-revolution Earth-orbit transfers. The parameterized control law in the form of the true optimal control is adopted and the control-related parameters are interpolated through a number of nodes along the time axis. The optimal control problem is thus converted into a parameter optimization problem that is in turn solved with nonlinear programming. Multiple shooting of mean elements and the technique for providing initial guesses through the Lyapunov control are employed to improve the convergence robustness of nonlinear programming. The parameter optimization problem thus formed could easily include any other static system parameters and possible equality and inequality constraints. The proposed direct optimization method serves as an ideal tool for solving low-thrust many-revolution transfer problems.

References

[1] Edelbaum T N. Propulsion requirements for controllable satellites. American Rocket Society Journal 1961; 31(8): 1079-1089.

[2] Sackett L L, Malchow H L, Edelbaum T N. Solar electric geocentric transfer with attitude constraints: analysis. NASA CR-134927, 1975.

[3] Betts J T. Very low-thrust trajectory optimization using a direct SQP method. Journal of Computational & Applied Mathematics 2000; 120(1): 27-40.

[4] 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.

[5] 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.

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

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

[8] Ilgen M R. Low thrust OTV guidance using Lyapunov optimal feedback control techniques. AAS/AIAA As-trodynamics Specialist Conference, Paper AAS 93-680. 1993; 1527-1545.

[9] Chang D E, Chichka D F, Marsden J E. Lyapunov functions for elliptic orbit transfer. Discrete and Continuous Dynamical System-Series B 2002; 2(1): 57-67.

[10] Petropoulos A E. Low-thrust orbit transfers using candidate Lyapunov functions with a mechanism for coasting. AIAA-2004-5089, 2004.

[11] Haberkorn T, Martinon P, Gergaud J. Low thrust minimum-fuel orbital transfer: a homotopic approach. Journal of Guidance, Control, and Dynamics 2004; 27(6): 1046-1060.

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

[13] Cefola P J, Long A C, Holloway G, Jr. The long-term prediction of artificial satellite orbits. AIAA-74-170, 1974.

[14] Danielson D A, Sagovac C P, Neta B, et al. Semiana-lytic satellite theory. NPS-MA-95-002, 1995.

[15] Lewis F L. Optimal control. New York: John Wiley & Sons, 1986; 149-154.

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

[17] Betts J T. Survey of numerical methods for trajectory optimization. Journal of Guidance, Control, and Dynamics 1998; 21(2): 193-207.

[18] Standish E M. The JPL planetary and lunary ephemerides DE402/LE402. Bulletin of the American Astronomical Society 1995; 27: 1203.

Biography:

Gao Yang Born in 1974, he received B.S. and M.S. degrees from Beijing University of Aeronautics and Astronautics and Chinese Academy of Sciences in 1997 and 2000, respectively. He received his Ph.D. degree in 2003 from University of Missouri-Columbia, USA. From 2004 to 2005, he was a postdoctoral research fellow in University of Missouri- Columbia. Since June 2005, he has been an associate research fellow at Academy of Opto-Electronic, Chinese Academy of Sciences. His academic interests include orbital mechanics, spacecraft dynamics and control, guidance and navigation, trajectory optimization, and space mission design.

E-mail: gaoy@aoe.ac.cn