A Numerical Study of Low-Thrust Limited Power Trajectories between Coplanar Circular Orbits in an Inverse-Square Force Field

A numerical study of optimal low-thrust limited power trajectories for simple transfer no rendezvous between circular coplanar orbits in an inverse-square force field is performed by two different classes of algorithms in optimization of trajectories. This study is carried out by means of a direct method based on gradient techniques and by an indirect method based on the second variation theory. The direct approach of the trajectory optimization problem combines the main positive characteristics of two well-known direct methods in optimization of trajectories: the steepest-descent first-order gradient method and a direct second variation second-order gradient method. On the other hand, the indirect approach of the trajectory optimization problem involves two different algorithms of the well-known neighboring extremals method. Several radius ratios and transfer durations are considered, and the fuel consumption is taken as the performance criterion. For small-amplitude transfers, the results are compared to the ones provided by a linear analytical theory.


Introduction
The main purpose of this work is to present a numerical study of optimal low-thrust limited power trajectories for simple transfers no rendezvous between circular coplanar orbits in an inverse-square force field.This study has been motivated by the renewed interest in the use of low-thrust propulsion systems in space missions verified in the last two decades due to the development and the successes of space missions powered by ionic propulsion; for instance, Mathematical Problems in Engineering Deep Space One and SMART 1 missions.Several researchers have obtained numerical and sometimes analytical solutions for a number of specific initial orbits and specific thrust profiles 1-6 .Averaging methods are also used in such researches 7-11 .Two idealized propulsion models have most frequently been used in the analysis of optimal space trajectories 12 : the limited power variable ejection velocity systems-LP systems-are characterized by a constraint concerning with the power there exists an upper constant limit for the power , and the constant ejection velocity-limited thrust systems-CEV systems-are characterized by a constraint concerning with the magnitude of the thrust acceleration which is bounded.In both cases, it is usually assumed that the thrust direction is unconstrained.The utility of these idealized models is that the results obtained from them provide good insight into more realistic problems.
In the study presented in this paper only LP systems are considered.The fuel consumption is taken as the performance criterion and it is calculated for various radius ratios ρ r f /r 0 , where r 0 is the radius of the initial circular orbit O 0 , and r f is the radius of the final circular orbit O f and for various time of flight t f − t 0 .The optimization problem associated to the space transfer problem is formulated as a Mayer problem of optimal control with Cartesian elements-components of position and velocity vectors-as state variables.Transfers with small, moderate, and large-amplitudes are studied, and the numerical results are compared to the results provided by a linear theory given in terms of orbital elements 12-15 .Two different classes of algorithms are applied in determining the optimal trajectories.They are computed through a direct approach of the trajectory optimization problem based on gradient techniques, and through an indirect approach based on the solution of the twopoint boundary value problem obtained from the set of necessary conditions for optimality.
The direct approach involves a gradient-based algorithm which combines the main positive characteristics of the steepest-descent first-order gradient method and of a direct method based upon the second variation theory second-order gradient method .This algorithm has two distinct phases: in the first one, it uses a simplified version of the steepestdescent method developed for a Mayer problem of optimal control with free final state and fixed terminal times, in order to get great improvements of the performance index in the first iterates with satisfactory accuracy.In the second phase, the algorithm switches to a direct method based upon the second variation theory developed for a Bolza problem with fixed terminal times and constrained initial and final states, in order to improve the convergence as the optimal solution is approached.This kind of algorithm for determining optimal trajectories is well known in the literature 16 , and the version used in this paper is quite simple, since it uses a simplified version of the steepest-descent method, as mentioned before, with terminal constraints added to the performance index by using a penalty function method see Section 2.2 .This procedure simplifies the algorithm, providing a solution with satisfactory accuracy, and can avoid some of typical divergence troubles of the classical steepest-descent method as discussed in McDermott and Fowler 17 .
The indirect approach involves the solution of the two-point boundary value problem through two different algorithms of the neighboring extremals method.The formulation of the neighboring extremals method, as presented herein, is associated with a Bolza optimal control problem with fixed initial and final times, fixed initial state and constrained final state 18, 19 .Basically, the method consists in iteratively determining the initial values of the adjoint variables and the Lagrange multipliers associated to the final constraints.It involves the linearization, about an extremal solution, of the nonlinear two-point boundary value problem defined by the application of Pontryagin Maximum Principle 20 to the optimization problem.The linearized problem has been solved through the state transition matrix, and through the generalized Riccati transformation 16, 27 .The algorithms based on the state transition matrix and the Riccati transformation are well known in the literature, and the version used in this paper has a slight modification as described in da Silva Fernandes and Golfetto 15 .
A brief description of the versions of the algorithms used in this paper can be found in da Silva Fernandes 21 .Finally, note that the results presented herein complete and extend the results previously obtained 15, 21, 22 .

Optimal Low-Thrust Limited Power Trajectories
In this section, the optimization problem concerning with optimal low-thrust limited power trajectories is formulated.Application of each one of the proposed algorithms is also presented.For completeness, a very brief description of the linear theory is included.

Formulation of the Optimization Problem
Low-thrust limited power propulsion systems are characterized by low-thrust acceleration level and a high specific impulse 12 .The ratio between the maximum thrust acceleration and the gravity acceleration on the ground, γ max /g 0 , is between 10 −4 and 10 −2 .For such system, the fuel consumption is described by the variable J defined as where γ is the magnitude of the thrust acceleration vector γ , used as control variable.The consumption variable J is a monotonic decreasing function of the instantaneous mass m of the space vehicle: where P max is the maximum power, and m 0 is the initial mass.The minimization of the final value J f is equivalent to the maximization of m f or the minimization of the fuel consumption.
The optimization problem concerning with simple transfers no rendezvous between coplanar orbits is formulated as: at time t, the state of a space vehicle M is defined by the radial distance r from the center of attraction, the radial and circumferential components of the velocity, u and v, and the fuel consumption J.In the two-dimensional formulation, the state equations are given by 23 : where μ is the gravitational parameter, R and S are the radial and circumferential components of the thrust acceleration vector, respectively.The optimization problem is stated as: it is proposed to transfer a space vehicle M from the initial state at the time t 0 0: to the final state at the prescribed final time t f : such that J f is a minimum, that is, the performance index is defined by: IP J t f .

2.6
Equations 2.4 and 2.5 are given in canonical units, and they define the initial and final circular orbits.u t f , v t f , and r t f denote the state variables at the prescribed final time t f , and 0, μ/r f , and r f are the prescribed values defining the final circular orbit.Similar definition applies at the initial time t 0 0 see 2.4 .For LP system, it is assumed that there are no constraints on the thrust acceleration vector 12 .
In the formulation of the optimization problem described above, the variables are written in canonical units, such that the gravitational parameter μ is equal to 1.

Application of the Gradient-Based Algorithm
As described in da Silva Fernandes 21 , the first phase of the gradient-based algorithm involves a simplified version of the steepest-descent method, which has been developed for a Mayer problem of optimal control with free final state and fixed terminal times.So, the optimal control problem defined by 2.3 -2.6 must be transformed into a new optimization problem with final state completely free.In order to do this, the penalty function method 24, 25 is applied.The new optimal control problem is then defined by 2.3 and 2.4 , with the new performance index obtained from 2.5 and 2.6 : where k 1 , k 2 , k 3 1.The penalty function method involves the progressive increase of the penalty constants; but, for simplicity, they are taken fixed in the gradient-based algorithm, since the steepest-descent is used only to provide a convex nominal solution as starting solution for the second order gradient method.
According to the algorithm of the simplified version of the steepest-descent method, the adjoint variables λ u , λ v , λ r , and λ J are introduced, and the Hamiltonian H is formed using 2. 3 21, 26 : From the Hamiltonian H, one finds the adjoint equations: and, from the performance index defined by 2.7 , we get the terminal conditions for the adjoint equations:

2.10
The algorithm also requires the partial derivatives of the Hamiltonian H with respect to the control variables.These partial derivatives are given by: The second phase of the gradient-based algorithm involves the second order gradient method, developed for a Bolza problem with fixed terminal times and constrained initial and final states, which requires the computation of the first order derivatives of the vector function ψ containing the terminal constraints and the scalar function Φ, corresponding to the augmented performance index, and the computation of the second order derivatives of the Hamiltonian H with respect to all arguments.The partial derivatives of the Hamiltonian function are given in a matrix form by: 12 α denotes the control vector α T R S , and it has been introduced to avoid confusion with the state variable u, x is the state vector x T u v r J , and λ is the adjoint vector λ T λ u λ v λ r λ J .From 2.5 and 2.6 , one defines the functions ψ and Φ: 13 where Λ i , i 1, 2, 3 are Lagrangian multipliers associated to the final constraints defined by 2.13 .The partial derivatives of ψ and Φ are then given by:

2.15
The results of the gradient-based algorithm to the optimization problem described above are presented in Section 3.

Application of the Neighboring Extremals Algorithms
Let us to consider the Hamiltonian function defined by 2.8 .Following the Pontryagin Maximum Principle 20 , the control variables R and S must select from the admissible controls such that the Hamiltonian function reaches its maximum along the optimal trajectory.Thus, The adjoint variables λ u , λ v , λ r , and λ J must satisfy the adjoint differential equations and the transversality conditions.Therefore, from 2.3 -2.5 and 2.16 one finds the following two-point boundary value problem for the transfer problem defined by 2.3 -2.6 : with the boundary conditions:

2.18
The neighboring extremals algorithms are based on the solution of a linearized twopoint boundary value problem that involves the derivatives of the right-hand side of 2.3 with respect to the state and adjoint variables 16, 21, 27 .These equations can be put in the following form: with δx t x n 1 t − x n t and δλ t λ n 1 t − λ n t , where n denotes the iterate, and A, B, and C are matrices given by: where

2.21
The results of the neighboring extremals algorithms to the optimization problem described above are presented in Section 3.

Linear Theory
For completeness, a very brief description of a first-order analytical solution for the problem of optimal simple transfer no rendezvous between close quasicircular coplanar orbits in an inverse-square force field is presented.This approximate solution, also referred as linear theory, is expressed in nonsingular orbit elements, and it is valid for orbits with very small eccentricities.According to Marec 12 or da Silva Fernandes and Golfetto 15 , for transfers between circular orbits, only Δα is imposed, and assuming that the initial and final positions of the vehicle in orbit are symmetric with respect to x-axis of the inertial reference system, that is, f − 0 Δ /2, the linear solution can be written as: where α a/a, h e cos ω, k e sin ω, where a is the semimajor axis, e is the eccentricity, ω is the argument of the pericenter, 0 n t − t 0 , and n μ/a 3 is the mean motion.The overbar denotes the reference orbit O about which the linearization is done.
The optimal thrust acceleration Γ * during the maneuver is expressed by: where e r and e s are unit vectors extending along radial and circumferential directions in a moving reference frame, respectively.The linear theory is applicable only for orbits which are not separated by large radial distance.If the reference orbit is chosen in the conventional way, that is, with the semimajor axis as the radius of the initial orbit, the radial excursion to the final orbit will be maximized 14 .A better reference orbit is defined with a semimajor axis given by an intermediate value between the values of semimajor axes of the terminal orbits.In this study, a is taken as a a 0 a f /2 in order to improve the accuracy in the calculations.
In the next section, the results of this linear theory are compared to the ones provided by the proposed algorithms.

Results
The results of a numerical analysis for optimal low-thrust limited power simple transfers no rendezvous between coplanar circular orbits in an inverse-square force field, obtained through the analytical and numerical methods described in the preceding sections, are presented for various radius ratios ρ r f /r 0 and for various time of flight t f − t 0 presented in Tables 1-8.All results are presented in canonical units as described in Section 2. A preliminary analysis of some interplanetary missions considering transfers from Earth to Venus, Mars, asteroid belt, Jupiter, and Saturn, which correspond to ρ 0.727, 1.523, 2.500, 5.203, and 9.519, respectively, is presented.In this preliminary analysis of interplanetary missions, the following assumptions are considered: 1 the orbits of the planets are circular; 2 the orbits of the planets lie in the plane of the ecliptic; 3 the flight of the space vehicle takes place in the plane of the ecliptic; 4 only the heliocentric phase is considered, that is, the attraction of planets on the spacecraft is neglected.show the values of the consumption variable J for small-amplitude transfers computed through the different approaches and the absolute relative difference in percent between the numerical and analytical results, according to the following definition:  The results provided by the neighboring extremals algorithm based on the state transition matrix denoted by number 1 have been chosen as the exact solution for each maneuver, in view of the accuracy obtained in fulfillment of the terminal constraints.Similar results for interplanetary transfers are presented in Tables 5, 6 and 7. Results for large-amplitude transfers with long time of flight are presented in Table 8.In both cases, the transfers are only computed through the neighboring extremals algorithms.

3.1
From the results presented in Tables 1-8, major comments are as follows: 1 The linear theory provides a very good approximation for the fuel consumption considering small-amplitude transfers with |ρ − 1| ≤ 0.100, that is, for transfers between close circular coplanar orbits.For the most of the maneuvers, d rel 1 < 0.5%; 2 For transfers with small time of flight t f − t 0 2.0, 3.0, 4.0, 5.0 time units , Tables 1, 2, 5 and 6 show that the maximum absolute relative difference d rel 1 occur for the most of the transfers with t f − t 0 5.This maximum value of d rel 1 is about 2% for ρ > 1 and 5.5% for ρ < 1;  5 For transfers between close orbits with small time of flight, the fuel consumption can be greatly reduced if the duration of the transfer increases: for instance, the fuel consumption for transfers with time of flight t f − t 0 2.0 time units is approximately ten times the fuel consumption for transfers with time of flight t f − t 0 4.0 time units , and, it is approximately hundred times the fuel consumption for transfers with time of flight t f − t 0 20.0 time units , considering any value of ρ; 6 For transfers with moderate time of flight, the fuel consumption decreases almost linearly with the time of flight; 7 For transfers with moderate amplitude ρ 0.727 and ρ 1.523 , Tables 5 and 6 show that 7.0% > d rel 1 > 1.0%; 8 Tables 1-6 show that the results obtained through the numerical algorithmsgradient and neighboring extremals-are very quite similar, regardless the amplitude of maneuver and the time of flight; 9 Table 7 shows that an Earth-asteroid belt mission with t f − t 0 20.0 time units approximately, 3.2 years and an Earth-Jupiter mission with t f − t 0 50.0 time units approximately, 8.0 years spend almost the same quantity of fuel.This result is closely related to the concept of transversals payoff curves in a field of extremals introduced by Edelbaum 7 in the study of optimal limited-power transfers in strong gravity field.Low-thrust limited power transfers with different amplitude ρ and different time of flight t f − t 0 can be performed with the same amount of fuel as shown in Figures 7, 8, and 9. From Figure 7, one finds that an Earth-Venus mission with t f − t 0 30.0 time units approximately, 4.8 years and an Earth-Venus mission with t f − t 0 36.5 time units approximately, 5.8 years also spend almost the same quantity of fuel, J 4.98 × 10 −4 canonical units .
In order to follow the evolution of the optimal thrust acceleration vector during the maneuver, it is convenient to plot the locus of its tip in the moving frame of reference.Figures 1, 2, 3, and 4 illustrate these plots for ρ 0.727, 0.950, 0.975, 1.025, 1.050, and 1.523, with t f − t 0 2.0, 3.0, 30.0 and 50.0.Note that the agreement between the numerical and analytical     results is better for small-amplitude transfers.For moderate-and large-amplitude transfers, this difference increases with the duration of the maneuvers.Similarly, Figures 5 and 6 show the evolution of the optimal thrust acceleration vector for large-amplitude transfers with long time of flight.Note that the magnitude of the thrust acceleration becomes smaller as the time of flight increases.
Figures 3, 4, and 5 show that for transfers with long time of flight the circumferential thrust acceleration is the main component of the optimal thrust acceleration.As the time of   flight increases the contribution of the radial component of the optimal thrust acceleration decreases and the optimal thrust acceleration tends to the circumferential acceleration.
Figures 7 and 8 show the consumption variable for interplanetary transfers with moderate flight of time, and Figure 9 shows the consumption variable for large-amplitude transfers with long time of flight.Note that transfers with different amplitude ρ and different time of flight t f − t 0 can be performed with the same amount of fuel see Comment 9 .

Conclusion
In this paper, a gradient-based algorithm and two different algorithms of the neighboring extremals method are applied to the analysis of optimal low-thrust limited power transfers between circular coplanar orbits in an inverse-square force field.The numerical results given by these algorithms have been compared to the analytical results obtained by a linear theory.The good agreement between these results shows that the linear theory provides a good approximation for the solution of the transfer problem concerned with small-amplitudes, that is, for transfers between close circular coplanar orbits.The linear theory can be used in preliminary mission analysis involving such kind of transfers.The results presented in the paper also show that the fuel consumption can be reduced if the time of flight of the transfer increases.For transfers with long time of flight, the circumferential thrust acceleration becomes the main component of the optimal thrust acceleration.A preliminary analysis of some interplanetary missions is presented using the neighboring extremals algorithms.

Figure 7 :Figure 8 :
Figure 7: Consumption for Earth-Venus and Earth-Mars transfers with moderate time of flight.

Table 1 :
Consumption variable J ρ > 1 for transfers with small time of flight.

Table 2 :
Consumption variable J ρ < 1 for transfers with small time of flight.

Table 3 :
Consumption variable J ρ > 1 for transfers with moderate time of flight.

Table 4 :
Consumption variable J ρ < 1 for transfers with moderate time of flight.

Table 5 :
Consumption variable J for Earth-Venus transfers.

Table 6 :
Consumption variable J for Earth-Mars transfers.

Table 7 :
Consumption variable J for interplanetary transfers with large-amplitude.

Table 8 :
Consumption variable J for large transfers.