Numerical and Analytical Study of Optimal Low-Thrust Limited-Power Transfers between Close Circular Coplanar Orbits

A numerical and analytical study of optimal low-thrust limited-power trajectories for simple transfer (no rendezvous) between close circular coplanar orbits in an inversesquare force field is presented. The numerical study is carried out by means of an indirect approach of the optimization problem in which the two-point boundary value problem, obtained from the set of necessary conditions describing the optimal solutions, is solved through a neighboring extremal algorithm based on the solution of the linearized twopoint boundary value problem through Riccati transformation. The analytical study is provided by a linear theory which is expressed in terms of nonsingular elements and is determined through the canonical transformation theory. The fuel consumption is taken as the performance criterion and the analysis is carried out considering various radius ratios and transfer durations. The results are compared to the ones provided by a numerical method based on gradient techniques.


Introduction
The main purpose of this paper is to present a numerical and analytical study of optimal low-thrust limited-power trajectories for simple transfers (no rendezvous) between close circular coplanar orbits in an inverse-square force field.The study of these transfers is particularly interesting because the orbits found in practice often have a small eccentricity and the problem of slight modifications (corrections) of these orbits is frequently met [1].Besides, the analysis has been motivated by the renewed interest in the use of low-thrust propulsion systems in space missions verified in the last two decades.Several researchers have obtained numerical, and sometimes analytical, solutions for a number of specific initial orbits and specific thrust profiles [2][3][4][5][6][7][8][9][10].Averaging methods are also used in such researches [11][12][13][14][15].
Low-thrust electric propulsion systems are characterized by high specific impulse and low-thrust capability and have a great interest for high-energy planetary missions and certain Earth orbit missions.For trajectory calculations, two idealized propulsion models are of most frequent use [1]: LP and CEV systems.In the power-limited variable ejection velocity systems or, simply, LP systems, the only constraint concerns the power, that is, there exists an upper constant limit for the power.In the constant ejection velocity limited-thrust systems or, simply, CEV systems, the magnitude of the thrust acceleration 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 this paper, only LP systems are considered.
In the study presented in the paper, 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 transfer durations t f − t 0 .Transfers with small and moderate amplitudes are considered.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.
The numerical study is carried out by a neighboring extremal algorithm which is based on the linearization about an extremal solution of the nonlinear two-point boundary value problem defined by the set of necessary conditions for a Bolza problem of optimal control with fixed initial and final times, fixed initial state, and constrained final state [16,17].The resulting linear two-point boundary value problem is solved through Riccati transformation.As briefly described in Section 2, a slight modification is introduced in the algorithm to improve the convergence.On the other hand, the analytical study is based on a linear theory expressed in terms of nonsingular orbital elements, similar to the ones presented in [1,18].Here, the linear theory is determined through canonical transformation theory using the concept of generalized canonical systems.This approach provides a simple way to compare the numerical solutions and the analytical theory.The numerical and analytical results are compared to the ones obtained through an algorithm based on gradient techniques [19,20].

Neighboring extremal algorithm based on Riccati transformation
For completeness, a brief description of the neighboring extremal algorithm used in the paper is presented in this section.This algorithm has a slight modification when compared to the well-known algorithms in the literature [16,17,21]: a constraint on the control variations is introduced.Numerical experiments have shown that this simple device improves the convergence.
Let the system of differential equations be defined by where x is an n-vector of state variables and u is an m-vector of control variables.It is assumed that there exist no constraints on the state or control variables.The problem consists in determining the control u * (t) that transfers the system (2.1) from the initial conditions to the final conditions at t f , and minimizes the performance index and ψ(•) : R n → R q , q < n, are assumed to be twice continuously differentiable with respect to their arguments.Furthermore, it is assumed that the matrix [∂ψ/∂x] has a maximum rank.By applying the Pontryagin maximum principle [21,22] to the Bolza problem with constrained final state and fixed terminal times defined by (2.1)-(2.4), the following twopoint boundary value problem is obtained: ) ) with ) ) where H(x,λ,u) = −F(x, u) + λ T f (x,u) is the Hamiltonian function, λ is an n-vector of adjoint variables and μ is a q-vector of Lagrange multipliers.The quantities H x ,H u ,g x ,..., and so forth, denote the partial derivatives.If x, λ, and u are taken to be column vectors, then H x , H λ , and H u are row vectors.In this way, ψ x is a q × n-matrix.The superscript T denotes the transpose of a matrix or a (row or column) vector.
Neighboring extremal methods are iterative procedures used for solving the two-point boundary value problem defined through (2.5)-(2.10).These methods are based on the second variation theory and consist in determining iteratively the unknown adjoint variables λ(t 0 ) and Lagrange multipliers μ.Let λ 0 (t 0 ) be an arbitrary starting approximation of the unknown adjoint variables at t 0 .The trajectory x 0 (t) corresponding to these starting values is obtained by integrating (2.5) from t 0 to t f , with the initial conditions (2.8).The vector of Lagrange multipliers μ is then calculated such that the transversality condition (2.9) is fulfilled.Since ψ x has a maximum rank, one finds that (2.11) Let λ 1 (t 0 ) = λ 0 (t 0 ) + δλ(t 0 ) and μ 1 = μ 0 + δμ be the next approximation.Following [16,21], the corrections (perturbations) δλ(t 0 ) and δμ are obtained in order to satisfy the linear two-point boundary value problem obtained from the linearization of (2. ) ) where the constant k, 0 < k ≤  [16,17,21].This accessory minimum problem is obtained expanding the augmented performance index, which includes through adjoint variables the constraints represented by the state equations, to second order and all constraints to first order, as described in the appendix.
According to the appendix, (2.14) must be replaced by (A.8) since a constraint on the control variations is imposed [23].Assuming that W 2 is chosen such that H uu + W 2 is nonsingular for t ∈ t 0 ,t f , we may solve (A.8) for δu(t), in terms of δx(t) and δλ(t): Substituting this equation into (2.12) and (2.13), it follows that ) where matrices A, B, and C are given by (2.21) We recall that the matrices A, B and C are evaluated on a nominal extremal solution.Equations (2.15) through (2.20) represent a linear two-point boundary value problem, whose solution can be obtained through a backward sweep method which uses the Riccati transformation [21]: where R is an n × n symmetric matrix, L is an n × q matrix, and Q is a q × q symmetric matrix.For (2.22) to be consistent with (2.15)-(2.20), the Riccati coefficients must satisfy the differential equations (2.23) with the boundary conditions (2.24) defined below.
The step-by-step computing procedure to be used in the neighboring extremal algorithm is summarized as follows.

Optimal low-thrust trajectories
In what follows, the neighboring extremal algorithm presented in previous section is applied to determine optimal low-thrust limited-power transfers between close coplanar circular orbits in an inverse-square force field.

Problem formulation.
A low-thrust limited-power propulsion system, or LP system, is characterized by low-thrust acceleration level and high specific impulse [1].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 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 of the fuel consumption J f is equivalent to the maximization of m f .The optimization problem concerned with simple transfers (no rendezvous) between coplanar orbits will be formulated as a Mayer problem of optimal control by using Cartesian elements as state variables.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 transverse components of the velocity, u and v, and the fuel consumption J. (Note that the radial component u should not be confused with the control variables defined in Section 2.) The geometry of the transfer problem is illustrated in Figure 3.1.
In the two-dimension optimization problem, the state equations are given by S. da Silva Fernandes and W. A. Golfetto 7 where μ is the gravitational parameter (it should not be confused with Lagrange multiplier defined in Section 2), R and S are the components of the thrust acceleration vector in a moving frame of reference, that is, Γ = Re r + Se s , with the unit vector e r pointing radially outward and the unit vector e s perpendicular to e r in the direction of the motion and in the plane of orbit.The optimization problem is stated as follows: it is proposed to transfer a space vehicle M from the initial conditions at t 0 , to the final state at the prescribed final time t f , such that J f is a minimum.We note that in the formulation of the boundary conditions above all variables are dimensionless.In this case, μ = 1.

Two-point boundary value problem.
Following the Pontryagin maximum principle [21,22], the adjoint variables λ u , λ v , λ r , and λ J are introduced and the Hamiltonian function H(u,v,r,J,λ u ,λ v ,λ r ,λ J ,R,S) is formed using the right-hand side of (3.3): The control variables R and S must be selected from the admissible controls such that the Hamiltonian function reaches its maximum along the optimal trajectory.Thus, we find that (3.7) The variables λ u , λ v , λ r , and λ J must satisfy the adjoint differential equations and the transversality conditions (2.6) and (2.9).Therefore, from (3.3)-(3.7),we get the following two-point boundary value problem for the transfer problem defined by (3.3)-(3.5): with the boundary conditions given by (3.4) and (3.5), and the transversality condition (3.9)

Applying the neighboring extremal algorithm.
The matrices A, B, and C describing the linearized two-point boundary value problem in the neighboring extremal algorithm can be obtained straightforwardly from (3.8) by calculating the partial derivatives of the right-hand side of the equation with respect to the state and adjoint variables taking into account a diagonal weighting matrix W 2 as described in the next paragraph.These matrices are then given by S. da Silva Fernandes and W. A. Golfetto 9 where (3.11) In Section 5, we present the results obtained through the neighboring extremal algorithm for several ratios ρ = r f /r 0 , ρ = 0.727;0.800;0.900;0.950;0.975;1.025;1.050;1.100;1.200;1.523,and nondimensional transfer durations of 2, 3, 4, 5.We note that the Earth-Mars transfer corresponds to ρ = 1.523 and Earth-Venus to ρ = 0.727.The criterion adopted for convergence is a tolerance of 1.0 × 10 −8 in the computation of corrections (variations) of the initial value of the adjoint variables.In view of this convergence criterion, the terminal constraints are obtained with an error less than 1.0 × 10 −6 , which means that ψ(x(t f )) ≤ 1.0 × 10 −6 .All simulations consider the transfer from low orbit to high orbit, with starting approximation λ 0 (t 0 ) = (0.001 0.001 0.001 −1), attenuation factor k = 0.10 for ρ = 0.727 and 1.5236, and k = 0.15 for the other values of ρ, and a diagonal matrix W 2 such that W 211 = W 222 = −σ, with σ = 5.5 for all maneuvers, except ρ = 0.727 and 0.800, with t f − t 0 = 4 and 5.In these cases, σ = 2.5.

Linear theory
In this section, a first-order analytical solution for the problem of optimal simple transfer defined in Section 3.1 is presented.
The Hamiltonian function H * governing the extremal (optimal) trajectories can be obtained as follows.Since λ J is a first integral (3.8) and λ J (t f ) = −1, from the transversality condition (3.9), it follows that λ J (t) = −1.Thus, from (3.7), we find the optimal thrust acceleration Introducing these equations into (3.6), it results that In what follows, we consider the problem of determining an approximate solution of the system of differential equations governed by the Hamiltonian H * by means of the theory of canonical transformations.This analytical solution is obtained through the canonical transformation theory using the concept of generalized canonical systems [24,25].Now, consider the Hamiltonian function describing a null thrust arc in the two-dimensional formulation of the optimization problem defined in Section 3.1:

.3) 10 Mathematical Problems in Engineering
Note that H is obtained from (3.6) taking R = S = 0 and adding the last term concerning the differential equation of the angular variable θ, which defines the position of the space vehicle with respect to a reference axis in the plane of motion.This variable is important for rendezvous problems and plays no special role for simple transfer problems like the one considered here, but it is necessary to define the canonical transformations described below.In the transformation theory described in the next paragraphs, it is assumed that the Hamiltonian H * is augmented in order to include this last term, that is, Note that H is the undisturbed part of H * and plays a fundamental role in our theory.
According to the properties of generalized canonical systems, the general solution of the system of differential equations governed by the Hamiltonian H can be expressed in terms of a fast phase and is given by [24,25]: where p is the semi latus rectum, e is the eccentricity, ω is the pericenter argument, f is the true anomaly (fast phase), and (λ p ,λ e ,λ f ,λ ω ) are adjoint variables to (p,e, f ,ω).Equations The undisturbed Hamiltonian function H is invariant with respect to this canonical transformation.Thus, the undisturbed Hamiltonian is written in terms of the new variables as The general solution of the new differential equations governed by the new undisturbed Hamiltonian function H is closely related to the solution of the time flight equation in the two-body problem for elliptic, parabolic, and hyperbolic motions [25].For quasicircular motions, this solution is very simple, as described in the next paragraphs.
Equations (4.5) have singularities for circular orbits (e = 0).In order to avoid this drawback, a set of nonsingular elements is introduced.The relationships between the singular orbital elements and the nonsingular ones are given by These equations define a Lagrange point transformation and the Jacobian of the inverse of this transformation must be computed in order to get the relationships between the corresponding adjoint variables.Thus, we get  p,e, f ,ω,λ p ,λ e ,λ f ,λ (4.10) Substituting (4.8) and (4.9) into (4.5),we get (4.11) These equations are valid for all orbits and define a Mathieu transformation between the Cartesian elements and the nonsingular orbital elements.
For quasicircular orbits, with very small eccentricities, (4.11) can be greatly simplified if higher-order terms in eccentricity are neglected.Thus, where n = μ/a 3 is the mean motion and = ω + M is the mean latitude.We note that first-order terms in eccentricity are retained in the state variables in order to get a trajectory with better accuracy.For adjoint variables, this is unnecessary since λ a , λ h , and λ k are small quantities for transfers between close circular orbits, that is, for small amplitude transfers.
Introducing (4.12) into the expression for H * , we finally get For transfers between close circular coplanar orbits, an approximate solution of the system of differential equations governed by H * can be obtained through simple integrations if the system is linearized about a reference circular orbit O with semimajor axis a.This solution can be put in the form where Δx = [Δα Δh Δk] T denotes the imposed changes on nonsingular orbital elements (state variables), α = a/a, λ α = aλ a , λ 0 is the 3 × 1 vector of initial values of the adjoint variables, and A is a 3 × 3 symmetric matrix.The adjoint variables are constant and the matrix A is given by a αα a αh a αk a hα a hh a hk a kα a kh a kk where ) ) ) ) ) Subscript f stands for the final time, = 0 + n(t − t 0 ), and t 0 is the initial time.The linear solution, described by (4.14)-(4.21), is in agreement with the one presented in [1,18], where it is obtained through a different approach.
In view of (4.1) and (4.12), the optimal thrust acceleration Γ * is expressed by The variation of the consumption variable ΔJ during the maneuver can be obtained straightforwardly from (4.13) and (4.15) by integrating, from t 0 to t f , the differential equation (see (3.3), (4.1), and (4.13)) where H * γ denotes the part of the Hamiltonian H * concerned with the thrust acceleration.Thus, where a αα ,a αh ,...,a kk are given by (4.16)-(4.21);and λ α , λ h , and λ k are obtained from the solution of the linear algebraic system defined by (4.14).We recall that the extremal (optimal) trajectory is given by (4.12) with the nonsingular elements a, h, and k calculated from (4.14).
For transfers between circular orbits, only Δα is imposed.If it is assumed that the initial and final positions of the vehicle in orbit are symmetric with respect to x-axis of the inertial reference frame, that is, f = − 0 = Δ /2, the solution of the system (4.14) is given by We note that the linear theory is applicable only to orbits which are not separated by large radial distance, that is, to transfers between close orbits.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 [26].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 our analysis, we have chosen a = (a 0 + a f )/2 in order to improve the accuracy in the calculations.
In the next section, the results of the linear theory are compared to the ones provided by the neighboring extremal algorithm described in Sections 2 and 3.

Results
The values of the consumption variable J computed through the neighboring extremal algorithm and the ones provided by the linear theory and by a numerical method based on gradient techniques [19] are presented in Tables 5.1 and 5.2, and are plotted in Figures 5.1 and 5.2, as function of the radius ratio ρ of the terminal orbits for various transfer durations t = t f − t 0 .The absolute relative difference in percent between the numerical and analytical results is also presented in the tables according to the following definition:  low-thrust limited-power transfer between close circular coplanar orbits in an inversesquare force field.Figures 5.1 and 5.2 also show that the fuel consumption can be greatly reduced if the duration of the transfer is increased.The fuel consumption for transfers with duration t f − t 0 = 2 is approximately ten times the fuel consumption for a transfer with duration t f − t 0 = 4.
In order to follow the evolution of the optimal thrust acceleration vector during the transfer, it is convenient to plot the locus of its tip in the moving frame of reference.Figures 5.3 and 5.4 illustrate these plots for small amplitude transfers with ρ = 0.950 and 1.050, and for moderate amplitude transfers, Earth-Mars (ρ = 1.523), and Earth-Venus (ρ = 0.727) transfers, with t f − t 0 = 2.It should be noted that the agreement between the numerical and analytical results is better for small amplitude transfers.For moderate amplitude transfers, this difference increases with the duration of the maneuvers.
Figures 5.5 and 5.6 show the time history of the state variables-u, v, and r-for a small amplitude transfer, ρ = 1.050, and a moderate amplitude transfer, ρ = 1.523, with t f − t 0 = 2. Again, the agreement between the numerical and analytical results is better for small amplitude transfers.The results-consumption variable, thrust acceleration, and trajectory-provided by the gradient-based algorithm are quite similar to the ones provided by the neighboring extremal algorithm.
It should be noted that the numerical algorithms based on the second variation theory-gradient-based algorithm and the neighboring extremal algorithm-provide quite similar results.This fact leads us to suppose that the solutions provided by the both algorithms are really optimal in the sense of a local minimum for the consumption variable J, although the sufficiency conditions are not tested.Besides, we note that the Pontryagin maximum principle is a set of necessary and sufficient conditions for the linearized problem describing the transfers between close circular orbits [26].

Conclusion
In this paper, a numerical and analytical study of optimal low-thrust limited-power trajectories for simple transfer (no rendezvous) between close circular coplanar orbits in an inverse-square force field is presented.The numerical study is carried out by means of a neighboring extremal algorithm and the analytical one is based on linear theory obtained through canonical transformation theory, using the concept of generalized canonical systems.The numerical and analytical results have been compared to the ones obtained through a numerical method based on gradient techniques.The great agreement between the results shows that the linear theory provides a good approximation for the
.1) Tables 5.1 and 5.2 show that d rel1 < 2% for ρ > 1, and d rel1 < 5.5% for ρ < 1.The greater values corresponds to transfers with moderate amplitude ρ = 0.7270 and ρ = 1.5236.On the other hand, d rel2 < 0.04% for all cases.Tables 5.1 and 5.2, and Figures 5.1 and 5.2 show the good agreement between the results.Note that the linear theory provides a good approximation for the solution of the