Analytical Solution for Optimal Low-Thrust Limited-Power Transfers Between Non-Coplanar Coaxial Orbits

1. Departamento de Ciência e Tecnologia Aeroespacial Instituto Tecnológico de Aeronáutica – Departamento de Matemática – São José dos Campos/SP – Brazil Correspondence author: Francisco das Chagas Carvalho | Departamento de Ciência e Tecnologia Aeroespacial Instituto Tecnológico de Aeronáutica – Departamento de Matemática – Praça Marechal Eduardo Gomes, 50 | Vila das Acácias | CEP: 12228-900 | São José dos Campos/SP – Brazil | E-mail: fchagas.carvalho@gmail.com Received: Jan. 31, 2017 | Accepted: Jun. 01, 2017 Section Editor: Antonio Bertachini ABSTRACT: In this paper, an analytical solution for time-fixed optimal low-thrust limited-power transfers (no rendezvous) between elliptic coaxial non-coplanar orbits in an inverse-square force field is presented. Two particular classes of maneuvers are related to such transfers: maneuvers with change in the inclination of the orbital plane and maneuvers with change in the longitude of the ascending node. The optimization problem is formulated as a Mayer problem of optimal control with the state defined by semi-major axis, eccentricity, inclination or longitude of the ascending node, according to the class of maneuver considered, and a variable measuring the fuel consumption. After applying Pontryagin’s maximum principle and determining the maximum Hamiltonian, short periodic terms are eliminated through an infinitesimal canonical transformation. The new maximum Hamiltonian resulting from this canonical transformation describes the extremal trajectories for long duration transfers. Closed-form analytical solution is then obtained through Hamilton-Jacobi theory. For long duration maneuvers, the existence of conjugate points is investigated through the Jacobi condition. Simplified solution is determined for transfers between close orbits. The analytical solution is compared to the numerical solution obtained by integration of the canonical system of differential equations describing the extremal trajectories for some sets of initial conditions. Results show a great agreement between these solutions for the class of maneuvers considered in the analysis. The solution of the two-point boundary value problem of going from an initial orbit to a final orbit, based on the analytical solution, is also discussed.


INTRODUCTION
In the last thirty years, important space missions have made use of low-thrust propulsion systems.Th e two pioneer missions were NASA-JPL Deep Space 1 and ESA's SMART-1.Deep Space 1 was the fi rst interplanetary spacecraft to utilize Solar Electric Propulsion.It was developed by NASA in the New Millennium Program to testing new technologies for future Space and Earth Science Programs.It was launched on October 24, 1998, and its mission terminated on December 18, 2001, when its fuel supply exhausted.SMART-1 was the fi rst of a series of ESA's Small Missions for Advanced Research in Technology.It was used to test solar electric propulsion and other deep-space technologies.It was launched on September 27, 2003, and its mission ended on xx/xx 03/28 where: γ is the magnitude of the thrust acceleration vector γ, used as a control variable.The consumption variable J is a monotonic decreasing function of the mass m of the space vehicle, (1) (2) (3) (4) (5) (6) (7) where: P max is the maximum power and m 0 is the initial mass.The minimization of the final value of the fuel consumption, J f , is equivalent to the maximization of m f .
The motion of a space vehicle P, powered by a limited-power engine in an inverse-square force field, can be described by the well-known Gauss equations of Celestial Mechanics.These equations are given by Battin (1987): where: a, e, I, Ω, ω, M are the classical orbital elements: a is the semi-major axis, e is the eccentricity, I is the inclination of the orbital plane, Ω is the longitude of the ascending node, ω is the argument of periapsis, and M is the mean anomaly; and r is the radial distance.R, S and W are, respectively, radial, circumferential and normal component of the thrust acceleration.E is the eccentric anomaly, f is the true anomaly and n is the mean motion, n 2 a 3 = µ, with µ denoting the gravitational parameter.
According to Marec andVinh (1977, 1980), the optimal transfers between coaxial orbits are coaxial transfers, that is, the argument of periapsis does not change during the transfer, and two classes of maneuvers are related to the transfers between non-coplanar coaxial orbits: change in the inclination of the orbital plane and change in the longitude of the ascending node.For the first maneuver, ω = 0 o or ω = 180 o and Ω = 0 o , and, for the second one, ω = 90 o or 270 o and I = 90 o .A sketch of the geometric xx/xx 04/28 configuration of the initial orbit and the final orbit for these maneuvers is represented in Fig. 1 (only "visible" portion of the orbits above the reference plane is represented).Therefore, at time t, the state of a space vehicle P is defined by five variables: a, e, I or W, M and J.The inclination of the orbital plane or the longitude of the ascending node is used according to the class of maneuver considered in the analysis.
The optimization problem can be formulated as a Mayer problem of optimal control as follows: It is proposed to transfer the space vehicle from the initial state (a 0 , e 0 , I 0 or W 0 , M 0 , 0) at time t 0 to the final state (a f , e f , I f or W f , M f , J f ) at time t f , such that the final consumption variable J f is a minimum.The duration of the transfer t f -t 0 is specified.For simple transfers (no rendezvous) the mean anomaly at time t f is free.
So, taking into account the characterization of the two classes of maneuvers related to the transfers between non-coplanar coaxial orbits, described in the preceding paragraphs, one finds that the state equations are given by Gauss' equations for semi-major axis, eccentricity and mean anomaly, as defined by Eqs. 2, 3 and 7, respectively, and, by the following differential equations: The variable θ is introduced to denote the inclination of the orbital plane or the longitude of the ascending node according to the maneuver considered.The upper (lower) sign refers to the maneuver with change in the inclination of the orbital plane with ω = 0 o (ω = 180 o ) or maneuver with change in the longitude of the ascending node with ω = 90 o (270 o ).Note that Eq. 8 is derived straightforwardly form Eqs. 4 or 5 following the conditions described for each maneuver -maneuver with change in the inclination of the orbital plane or maneuver with change in the longitude of the ascending node -as described in a preceding paragraph.
According to the Pontryagin's maximum principle (Pontryagin et al. 1962), the adjoint variables (Lagrange multipliers associated to the constraints represented by state equations) p a , p e , p θ , p M and p J are introduced and the Hamiltonian function H (a, e, θ, M, J, p a , p e , p θ , p M , p J , R, S, W) is formed using the right-hand side of Eqs. 8 and 9, xx/xx 05/28 The control variables R, S and W must be selected from the admissible controls such that the Hamiltonian function reaches its maximum along the optimal trajectory.Thus, the components of the optimal thrust acceleration are given by: (10) These equations can be simplified if it will be taken into account that the adjoint variable p J is a first integral of the canonical system governed by the maximum Hamiltonian; that is, p J is a constant whose value is obtained from the transversality condition, p J (t J ) = -1.Accordingly, p J (t J ) = -1.So, introducing this result into Eqs.11 and 12, one finds the components of the optimal thrust acceleration.
Introducing Eqs.11-13 into the Eq. 10, one finds: where H 0 denotes the undisturbed Hamiltonian and H* γ is the part related to the optimal thrust acceleration.These functions are expressed by: It also be noted that the part of the maximum Hamiltonian related to the optimal thrust acceleration can be expressed in a compact form: This result is used to compute the fuel consumption by simple quadrature of Eq. 9.
The problem of determining a first-order analytical solution of the system of differential equations governed by the maximum Hamiltonian H* is solved by means of the theory of canonical transformations as it will be described later.

FIRST-ORDER ANALYTICAL SOLUTION
A first-order analytical solution of the canonical system governed by the maximum Hamiltonian H* is derived by means of canonical transformation theory.Firstly, the short periodic terms are eliminated through an infinitesimal canonical transformation built using Hori method (Hori 1966).Hori method is a perturbation method based on Lie commutation theorem (Lie 1888) and it provides an explicit canonical transformation between old and new canonical variables, differently from the classical perturbation methods based on Hamilton-Jacobi theory, in which the transformation of variables involves new and old variables in a mixed form and requires the solution of an inversion problem to get the final form of the transformation of variables.We note that the new maximum Hamiltonian resulting from the infinitesimal canonical transformation describes the extremal trajectories for long duration transfers.A set of first integrals is then derived for the canonical system of differential equations governed by the new Hamiltonian function and a closed-form analytical solution is obtained through Hamilton-Jacobi theory.The separation of variables technique is applied to solve the Hamilton-Jacobi equation (Lanczos 1970).

ELIMINATION OF SHORT PERIODIC TERMS
In order to eliminate the short periodic terms from the maximum Hamiltonian, it is assumed that the undisturbed Hamiltonian H 0 is of zero order and the disturbing part H * γ is of the first-order in a small parameter defined by the magnitude of the thrust acceleration.
Consider an infinitesimal canonical transformation, The new variables are designated by the prime.
According to the algorithm of Hori method, at order 0, one finds: F 0 denotes the new undisturbed Hamiltonian.In order to determine the higher order terms of the new Hamiltonian and the generating function of the canonical transformation, the algorithm proposed by da Silva Fernandes (2003) is applied.According to this algorithm, based on the method of variation of constants, the Poisson brackets involving the generating function and the new undisturbed Hamiltonian, that defines the general equation of the algorithm, is converted into a time partial differential equation.This procedure simplifies the determination of the new Hamiltonian and the generating function, since it eliminates the Hori auxiliary parameter.
So, consider the canonical system described by F 0 : xx/xx 07/28 general solution of which is given by: The subscript 0 denotes the constants of integration.Introducing this general solution into the equation of order 1 of the algorithm of Hori method, one finds: with the functions H* γ written in terms of the constants of integration.Following the proposed algorithm, one finds: Terms factored by p' M have been omitted in equations above, since only transfers (no rendezvous) are considered.Note that p' M is a first integral of the average canonical system; that is, p' M is constant and its value is computed from the transversality condition which gives p' M (t f )= 0 for simple transfers (the final position of the space vehicle is free).This result simplifies the expressions of the new Hamiltonian and the generating function as given by equations above.
The subscript 0 denotes the constants of integration.Introducing this general solution into the equation of order 1 of the algorithm of Ho method, one finds with the functions * H  written in terms of the constants of integration.
Following the proposed algorithm, one finds Terms factored by M p have been omitted in equations above, since only transfers ( rendezvous) are considered.Note that M p is a first integral of the average canonical system that is, M p is constant and its value is computed from the transversality condition which giv for simple transfers (the final position of the space vehicle is free).This resu xx/xx 08/28

A SET OF FIRST INTEGRALS OF THE AVERAGED CANONICAL SYSTEM
Consider the Mathieu transformation (Lanczos 1970) defined by the following equations: (25) The Hamiltonian function F 1 is invariant with respect to this transformation, The averaged canonical system governed by the Hamiltonian function F'' has a set of first integrals, which can be derived straightforwardly from the differential equations as described by Pines (1964) and Edelbaum and Pines (1970).
The first integrals are given by the following equations:

CLOSED-FORM SOLUTION OF THE CANONICAL SYSTEM GOVERNED BY HAMILTONIAN FUNCTION F''
Consider the canonical transformation, defined by a generating function W such that the constants C 1 , C 2 and E become the new generalized coordinates.
Since the new Hamiltonian function F'' is a quadratic form in the adjoint variables or conjugate momenta, the separation of variables technique will be applied for solving the Hamilton-Jacobi equation (Lanczos 1970).
Consider the transformation equations: xx/xx 09/28 where: is the generating function.
The corresponding Hamilton-Jacobi equation is then given by: Following the separation of variables technique (Lanczos 1970), it is assumed that the generating function W is equal to a sum of functions, each of which depends on a single old variable, that is: Therefore, from Eqs. 22-29, it follows that A complete solution of these equations is given by: with 5C 2 1 + C 2 2 = 5C 2 .We note that W 2 is given as indefinite integral, since only its partial derivatives are needed (see Eqs. 25-27), as shown in the next paragraphs.Now, consider the differential equations for the conjugate momenta of the canonical system governed by the new Hamiltonian F'' = E.These equations are: . .Introducing the generating function W, defined by Eqs.29-32, into the transformation equations, defined by Eqs.25-27, and taking into account the general solution defined by Eqs.33 for the conjugate momenta, one finds: From the above equations one finds, after some simplifications, the solution of the canonical system governed by the Hamiltonian F'' for a given set of initial conditions: This equation can also be used to evaluate the adjoint variable p'' a along the trajectory, as well as Eq.45.

xx/xx 12/28
An alternative expression for variable ψ can be determined from Eqs. 42, 49, 51, 52 and 53, as it follows: This equation will be useful in the analysis of conjugate point for long duration transfers, discussed later.Note that ψ does not depend on k 0 .

FIRST-ORDER SOLUTION OF THE CANONICAL SYSTEM GOVERNED BY HAMILTONIAN FUNCTION H*
Following the algorithm of Hori method (Hori 1966), the partial derivatives of the generating function S 1 with respect to the adjoint variables must be determined in order to obtain a first-order solution of the canonical system governed by Hamiltonian functionH*.Computing these partial derivatives and applying the initial conditions, one finds: with a' , e' , …, p' θ previously defined through the Mathieu transformation given by Eqs.19-21.The eccentric anomaly E' is computed from the well-known Kepler's equation with the mean anomaly given by: An approximate expression for the fuel consumption along the extremal trajectory is obtained by simple quadrature of Eq. 9 and can be put in a compact form as (da Silva Fernandes and das Chagas Carvalho 2008): with ΔS 1 = S 1 (M') -S 1 (M' 0 ), with S 1 given by Eq. 18.Note that Eq. 60 includes the periodic terms through the generating function S 1 .

LONG DURATION TRANSFERS
In what follows is presented a discussion about the existence of conjugate points through the analysis of Jacobi condition, and a description of an algorithm for solving the two-point boundary value problem of going form an initial orbit to a final orbit for long duration transfers.Such transfers are described by the general solution of the canonical system governed by the Hamiltonian F'' xx/xx 13/28 EXISTENCE OF CONJUGATE POINT Consider Eqs.41-44, which defi ne a two-parameter family of extremals in the phase space (a'' , φ'' , θ '') for a given initial phase point (a 0 '' , φ 0 , θ 0 '') corresponding to an initial orbit.By eliminating the auxiliary variables τ and ψ (see Eq. 55), α =a''/ a 0 '' and θ'' can be written as explicit functions of φ, that is, α = α(φ, φ 0 , k 0 , k 1 ) and θ'' (φ, φ 0 , θ 0 '' , k 1 ).Th e conjugate points to the phase point (a 0 '' , φ 0 , θ 0 '') are determined through the roots of the equation (Bliss 1946): Since θ'' does not depend on k 0 , ∂θ''/∂k 0 = 0. On the other hand, from Eq. 41 and from the remark about Eq. 55, one fi nds: Th us, the problem of determining the conjugate points reduces to the analysis of the roots of the following equation: From Eqs. 43 and 44, it follows that Th e partial derivative ∂θ'' = ∂k 1 is then given by Th e auxiliary variable τ is reintroduced to simplify the expression.
Th erefore, aft er some simplifi cations, one fi nds that the conjugate points are given by the roots of following equation: If a conjugate point τ* exists, or, equivalently, φ* (see Eq. 43 connecting τ and φ), then the extremal does not yield a local minimum.Th is is the necessary Jacobi condition in the Calculus of Variations for an extremum (Bliss 1946).

SOLUTION OF THE TWO-POINT BOUNDARY VALUE PROBLEM
In what follows, the solution of the two-point boundary value problem of going from an initial orbit to a fi nal orbit for long duration transfers is discussed.
By eliminating k 0 from the above equations, one finds: On the other hand, J = Et, thus: Note that t 0 = 0 in equations above.Equation 67 describes the time evolution of the fuel consumption during a long-time maneuver.
The solution of the two-point boundary value problem of going from an initial orbit O 0 : (a 0 , e 0 , θ 0 ) to a final orbit O f : (a f , e f , θ f ) defined by Eqs.41-44 involves the determination of the initial value of the adjoint variables p'' a0 , p'' φ0 and p'' θ0 , or equivalently, the determination of the auxiliary constants k 0 , k 1 and k 2 .Note that k 2 is written in terms of e 0 , θ'' 0 and k 1 through Eqs.43 and 50.Accordingly, the solution of the two-point boundary value problem reduces to determine the other two constants.
Note that from the preceding section, Eq. 61, θ'' only depends on k 1 .Accordingly, k 1 can be determined by solving Eq. 61, iteratively by means of Newton-Raphson method, for the given final conditions; that is, for θ''(t f ) = θ f .Thus, with the partial derivative (∂θ''/ ∂k 1 ) k 1 = k 1n computed from Eq. 62.Note that the Newton-Raphson algorithm fails, if a conjugate point exists.
The constant k 0 is then determined as follows.Given k 1 , the auxiliary variable ψ f can be calculated from Eq. 55.Solving Eq. 41 for k 0 , one finds: Once the two-point boundary value problem has been solved, the fuel consumption variable J can be calculated from Eq. 67.The initial value of the adjoint variables p'' a0 , p'' φ0 and p θ0 are then obtained as follows.From Eqs. 41 and 64, one finds: with u f /ν 0 given by Eq. 66.
Now, solving Eqs.48 and 49 for p'' φ0 and p'' θ0 , it follows that The steps of an iterative algorithm for solving the two-point boundary value problem can be described as:

•
For a given set of initial conditions (a 0 , e 0 ) and final conditions (a f , e f ), compute α f , φ 0 = sin -1 e 0 and φ f = sin -1 e f .

•
Guess a starting value for k 1 .

•
Compute successively p'' a0 , p'' φ0 and p'' θ0 using the Eqs.69-71.The solution of the two-point boundary value problem for long duration transfers is used as starting guess for the solution of the complete problem as described later.

TRANSFERS BETWEEN CLOSE ORBITS
For transfers between close non-coplanar coaxial orbits, a simplified and complete first-order solution can be obtained through a linearization of the right-hand side of Eqs.56-58, 59 and 60 about a reference orbit O with semi-major axis aand eccentricity For transfers between close orbits, the consumption variable can be written as: with a αα , a ee , and a θθ given by Eqs.74, 77 and 78, respectively.

SOLUTION OF THE TWO-POINT BOUNDARY VALUE PROBLEM
An iterative algorithm based on the complete first-order analytical solution is briefly described for solving the two-point boundary value problem of going from an initial orbit O 0 :(a 0 , e 0 , θ 0 ) to a final orbit O f :(a f , e f , θ f ) at the prescribed final time t f .
For a given final time, Eqs. 56-58, 59 and 60 can be represented as follows: xx/xx 17/28 where y 1 = a, y 2 = e and y 3 = θ.Note that p a0 , p e0 and p θ0 appear explicitly in the short periodic terms and also implicitly through a' (t), e' (t), θ'(t) and E'(t).Thus, functions g i , i = 1, 2, 3, are nonlinear in these variables.So, the two-point boundary value problem can be stated as: Find p a0 , p e0 and p θ0 such that y 1 (t f ) = a f , y 2 (t f ) = e f , y 3 (t f ) = θ f .This problem can be solved through a Newton-Raphson algorithm (Stoer and Bulirsch 2002), with the partial derivatives of functions g i computed numerically by means of a procedure of centered differences.As previously mentioned, the starting guess for the iterative procedure is given by the solution of the two-point boundary value problem concerning transfers with long duration.

RESULTS
Two types of problems are considered: a direct problem, which corresponds to generate extremals trajectories for a given set of initial conditions for the state and adjoint variables, and an indirect problem concerning the solution of the two-point boundary value problem.In the first problem, a comparison between the complete, nonlinear and linear, firstorder analytical solutions, derived in the preceding sections, and a numerical solution obtained by integrating the canonical system of differential equations governed by new Hamiltonian function H*, given by Eqs. 14 and 15, is discussed and the accuracy of the analytical approach is established.The second problem involves the comparison of the complete first-order analytical solution, including the short periodic terms, and the secular analytical solution (which describes the long duration transfers) in solving the two-point boundary value problem of going from an initial orbit O 0 :(a 0 , e 0 , θ 0 ) to a final orbit O f :(a f , e f , θ f ) at the prescribed final time t f .A Runge-Kutta-Fehlberg method of orders 4 and 5, with step-size control, relative error tolerance of 10 -11 , and absolute error tolerance of 10 -12 , as described in Forsythe et al. (1977), has been used in the numerical integrations.

DIRECT PROBLEM
Figures 2-4 show the results for a comparison between three distinct solutions: the complete first-order analytical solution, the secular analytical solution, and the numerical solution.Two sets of initial conditions (state and adjoint variables) and transfer duration, defined in Table 1, are used in the comparison.Note that a transfer with moderate time of flight and a long duration transfer are computed.In Table 2, final values of state variables are shown.In these tables and in all figures, the state variables are expressed in canonical units, except the inclination of the orbital plane given     in degrees.It should be noted that similar results can be obtained for maneuver with modification of the ascending node, since this maneuver is equivalent to the one with modification of the inclination of the orbital plane, taking into account the conditions described earlier (see paragraph after Eqs.2-7).For simplicity, the results refer to maneuvers with modification of the orbital plane.
From the results presented in Figs.2-4 and 5-7, and in Table 2, note that there exists an excellent agreement between the complete analytical solution and the numerical solution.The analytical solution has the same accuracy of the numerical solution for all state variables -semi-major axis, eccentricity, inclination of the orbital plane and consumption variable -considering moderate time of flight (t f -t 0 = 25.0) and large time of flight (t f -t 0 = 500.0).On the other hand, by analyzing the time evolution of the semi-major axis, eccentricity and inclination of the orbital plane represented in the plots of Figs.2-4 and 5-7, as given by the complete analytical solution, the numerical solution and the secular solution, we can see that the amplitudes of the short periodic terms are small, but they can be significant for transfers with short or moderate duration.Note that the difference between the numerical solution and the secular solution decrease as the transfer duration increases.The good agreement between the numerical solution and the complete analytical solution suggest that the latter can be xx/xx 20/28 used in the solution of the two-point boundary value problem of going from an initial orbit to a final orbit, as previously described.This last remark is relevant, since the numerical algorithm to compute optimal low-thrust trajectories based on the analytical solution requires a simple Newton-Raphson method to solve the two-point boundary value problem, as described above.

INDIRECT PROBLEM
In view of the previous results, one sees that the complete analytical first-order solution gives an accurate solution for the extremals trajectories concerning the transfers between coaxial non-coplanar orbits.So, in this section a study of some maneuvers is made using the analytical solution.The solution of the two-point boundary value problem of going from an initial orbit O 0 :(a 0 , e 0 , θ 0 ) to a final orbit xx/xx O f :(a f , e f , θ f ) at the prescribed final time t f is obtained considering the two versions of the analytical solution: the complete non-linear solution and the simplified linear solution for transfers between close orbits.This analysis allows us to discuss the applicability of the linear solution.Eight different maneuvers are considered, as described in Table 3.The initial orbit O 0 is the same for all maneuvers and it is defined by the following set of orbital elements: a 0 = 1.0 canonical unit, e 0 = 1.0 and I 0 = 1.0 degrees.The final value of the consumption variable J is taken as a comparison parameter for five transfer durations.Table 4 shows the results obtained through the two analytical solutions.J Non-linear is the consumption variable computed by Eq. 60, and J Linear is the consumption variable computed by Eq. 80.In both cases, variable J is expressed in canonical units.Figures 8-13   xx/xx the time evolution of the semi-major axis, the eccentricity, the inclination of the orbital plane and the consumption variable for maneuvers 1, 3 and 7. Two transfer durations are considered: t f -t 0 = 25.0 and t f -t 0 = 250.0time units.All variables are expressed in canonical units, except the inclination of the orbital plane, given in degrees.It should be noted that there exists a good agreement between the two analytical solutions, linear and non-linear, for maneuver 1, which corresponds to a transfer between close orbits.On the other hand, for maneuver 3, which corresponds to a transfer with moderate amplitude, the results given by the linear solution are not good.Similar conclusion applies to maneuver 7, although the behavior for larger duration transfer is fair.Now, consider the problem of solving the two-point boundary value problem for long-time transfers using the secular solution and the complete analytical solution, which includes the short periodic terms.In order to compare these solutions, maneuver 3 described in Table 3 is considered, with two transfer durations, t f -t 0 = 125.0 and 500.0time units.Table 5 shows the initial value of the adjoint variables obtained by using the secular solution and by using the complete solution.Note that there exist small differences between these initial values.As previously described, a small adjust of the initial value of the adjoint variables is made by a Newton-Raphson algorithm when the short periodic terms are included.Note that this adjust become smaller for a very large transfer duration.Moreover, show that the secular solution of the TPBVP does not represent a mean solution of the complete solution of the same TPBVP.
Figure 17 represents the optimal trajectory for maneuver 3 with time of flight equal to 125.0 time units.Note that the motion of the vehicle resembles a spiral, departing from the initial orbit and arriving at the terminal orbit after twelve revolutions around the central body.If the maneuver represents a transfer between two orbits around the Earth with the semi-major axis of the initial orbit, a 0 = 1.0 distance unit, corresponding to 8000 km, then the maneuver lasts 39.35 hours and the magnitude of the average acceleration is equal to 1.59 cm/s 2 .For such maneuver, the ratio between the average acceleration and the gravity acceleration on the ground is approximately 1.6 × 10 -3 , a typical value for low-thrust propulsion system.Similar conclusions are obtained for the other maneuvers discussed in this section.

CONCLUSION
An approximated first-order analytical solution for optimal time-fixed low-thrust limited power transfers (no rendezvous) between elliptic coaxial non-coplanar orbits in an inverse-square force field is determined through canonical transformation theory.The existence of conjugate point for long duration transfers is investigated through Jacobi condition.For transfers between close orbits a simplified solution is straightforwardly obtained by linearizing the new Hamiltonian and the generating function obtained through Hori method.For the direct problem of generating extremals trajectories, the analytical solution is compared to the numerical solution obtained by integrating the canonical system of differential equations describing the extremal trajectories for some sets of initial conditions, and the accuracy of the analytical approach is established.For the indirect problem, an iterative algorithm based on the first-order analytical solution is described for solving the two-point boundary value problem of going from an initial orbit to a final orbit.A comparison between the linearized and nonlinear analytical solutions is made and it has been noticed that the linearized solution can provide good results for transfers between close orbits.

Figure 1 .
Figure 1.Geometric configuration of the initial orbit and the final orbit: (a) represents the maneuver with change in the inclination of the orbital plane and (b) represents the maneuver with change in the longitude of ascending node.
with the auxiliary constants k 0 , k 1 and k 2 defined as functions of the initial value of the adjoint variables by: e -.This simplified solution can be put in the form:where Δx = [Δα Δe Δθ] T denotes the imposed changes on orbital elements (state variables), α = a/a -, p α = ap a , p 0 denotes the 3 × 1 vector of initial values of the adjoint variables, and, A is a 3 × 3 symmetric matrix.In this simplified solution, the adjoint xx/xx 16/28 variables are constant and the matrix A is given by: anomaly determined from Kepler's equation with the mean anomaly M -= M -0 + n -(t -t 0 ) and t 0 is the initial time.Since the imposed variations on the orbital elements are linear in the adjoint variables, Eq. 72, the solution of the two-point boundary value problem is very simple and can be obtained by standard techniques.

Figure 4 .
Figure 4. Comparison between secular, analytical and numerical time evolution of inclination of orbital plane for maneuver 1 with t f -t 0 = 25.0 and t f -t 0 = 500.0.

Figure 5 .
Figure5.Comparison between secular, analytical and numerical time evolution of semi-major axis for maneuver 2 with t f -t 0 = 25.0 and t f -t 0 = 500.0.

Figure 6 .
Figure 6.Comparison between secular, analytical and numerical time evolution of eccentricity for maneuver 2 with t f -t 0 = 25.0 and t f -t 0 = 500.0.

Figure 7 .
Figure 7.Comparison between secular, analytical and numerical time evolution of inclination of orbital plane for maneuver 2 with t f -t 0 = 25.0 and t f -t 0 = 500.0.

Figure 9 .
Figure 9.Time evolution of state variables for maneuver 1 with t f -t 0 = 250.0.

Figure 15 .
Figure 15.Secular solution and complete analytical solution of the TPBVP for maneuver 3 for both transfer durationsevolution of eccentricity.

Figure 16 .
Figure 16.Secular solution and complete analytical solution of the TPBVP for maneuver 3 for both transfer durationsevolution of inclination of orbital plane.

Table 1 .
Set of initial conditions and transfer duration.

Table 2 .
Final state variables. show

Table 3 .
Set of terminal orbits.
Time evolution of state variables for maneuver 1 with t f -t 0 = 25.0.