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

Authors: Sandro da Silva Fernandes; Carlos Roberto Silveira Filho; Wander Almodovar Golfetto
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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, 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ρ=rf/r0, where r0 is the radius of the initial circular orbit O0, and rf is the radius of the final circular orbit Of and for various time of flight tf-t0. 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 two-point 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 steepest-descent 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].


## 2. 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.


### 2.1. 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/g0, is between 10−4 and 10−2. For such system, the fuel consumption is described by the variable J defined as
(2.1)J=12∫t0tγ2dt,
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:
(2.2)J=Pmax⁡(1m-1m0),
where Pmax is the maximum power, and m0 is the initial mass. The minimization of the final value Jf is equivalent to the maximization of mf or the minimization of the fuel consumption.The optimization problem concerning with simple transfers (no rendezvous) between coplanar orbits is formulated as: at timet, 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]:(2.3)dudt=v2r-μr2+R,dvdt=-uvr+S,drdt=u,dJdt=12(R2+S2),
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 t0=0:
(2.4)u(0)=0v(0)=1r(0)=1J(0)=0,
to the final state at the prescribed final time tf:
(2.5)u(tf)=0v(tf)=μrfr(tf)=rf,
such that Jf is a minimum, that is, the performance index is defined by:(2.6)IP=J(tf).
Equations (2.4) and (2.5) are given in canonical units, and they define the initial and final circular orbits. u(tf),v(tf), and r(tf) denote the state variables at the prescribed final time tf, and 0, μ/rf, and rf are the prescribed values defining the final circular orbit. Similar definition applies at the initial time t0=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.


### 2.2. 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):(2.7)IP=J(tf)+k1(u(tf))2+k2(v(tf)-1rf)2+k3(r(tf)-rf)2,
where k1,k2,k3 ≫ 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]:(2.8)H=λu(v2r-μr2+R)+λv(-uvr+S)+λru+12λJ(R2+S2).From the HamiltonianH, one finds the adjoint equations:
(2.9)dλudt=vrλv-λr,dλvdt=-2vrλu+urλv,dλrdt=(v2r2-2μr3)λu-uvr2λv,dλJdt=0,
and, from the performance index defined by (2.7), we get the terminal conditions for the adjoint equations:
(2.10)λu(tf)=-2k1u(tf),λv(tf)=-2k2(v(tf)-1rf),λr(tf)=-2k3(r(tf)-rf),λJ(tf)=-1.The algorithm also requires the partial derivatives of the HamiltonianH with respect to the control variables. These partial derivatives are given by:(2.11)∂H∂R=λu+RλJ,∂H∂S=λv+SλJ.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:
(2.12)Hαα=[λJ00λJ],Hλα=[100100RS],Hxα=0(4×2nullmatrix),Hλx=[02vr-v2r2+2μr30-vr-uruvr2010000000],Hxx=[0-λvrvr2λv0-λvr2λur-2vr2λu+ur2λv0λvvr2-2vr2λu+ur2λv(2v2r3-6μr4)λu-2uvr3λv00000].α denotes the control vector αT=[RS], and it has been introduced to avoid confusion with the state variable u,x is the state vector xT=[uvrJ], and λ is the adjoint vector λT=[λuλvλrλJ].From (2.5) and (2.6), one defines the functions ψ  and Φ:(2.13)ψ=[u(tf)v(tf)-1rfr(tf)-rf],(2.14)Φ=J(tf)+Λ1u(tf)+Λ2(v(tf)-1rf)+Λ3(r(tf)-rf),
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)ψx=[100001000010],Φxx=0.(4×4nullmatrix).The results of the gradient-based algorithm to the optimization problem described above are presented in Section3.


### 2.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,
(2.16)R*=-λuλJ,S*=-λvλJ.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):
(2.17)dudt=v2r-μr2-λuλJ,dvdt=-uvr-λvλJ,drdt=u,dJdt=12λJ2(λu2+λv2),dλudt=vrλv-λr,dλvdt=-2vrλu+urλv,dλrdt=(v2r2-2μr3)λu-uvr2λv,dλJdt=0,
with the boundary conditions:
(2.18)u(0)=0,v(0)=1,r(0)=1,J(0)=0,u(tf)=0,v(tf)=μrf,r(tf)=rf,λJ(tf)=-1.The neighboring extremals algorithms are based on the solution of a linearized two-point boundary value problem that involves the derivatives of the right-hand side of (2.17) with respect to the state and adjoint variables [16, 21, 27]. These equations can be put in the following form:
(2.19)dδxdt=Aδx+Bδλ,dδλdt=Cδx-ATδλ,
with δx(t)=xn+1(t)-xn(t) and δλ(t)=λn+1(t)-λn(t), where n denotes the iterate, and A, B, and C are matrices given by:
(2.20)A=[02vr-a0-vr-uruvr2010000000],B=[-1λJ00λuλJ20-1λJ0λvλJ20000λuλJ2λvλJ20-cλJ3],C=[0λvr-vλvr20λvr-2λurb0-vλvr2bd00000],
where
(2.21)a=v2r2-2μr3,b=2vr2λu-ur2λv,c=λu2+λv2,d=(-2v2r3+6μr4)λu+2uvr3λv.The results of the neighboring extremals algorithms to the optimization problem described above are presented in Section3.


### 2.4. 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:
(2.22)u=na(hsin⁡l̅-kcos⁡l̅),v=na(1+hcos⁡l̅+ksin⁡l̅),r=a1+hcos⁡l̅+ksin⁡l̅,J(t)=12a̅5μ3{[52(l̅-l̅0)-34(sin⁡2l̅-sin⁡2l̅0)]λk24(l̅-l̅0)λα2+8(sin⁡l̅-sin⁡l̅0)λαλh-8(cos⁡l̅-cos⁡l̅0)λαλk+[52(l̅-l̅0)+34(sin⁡2l̅-sin⁡2l̅0)]λh2-32(cos⁡2l̅-cos⁡2l̅0)λhλk+[52(l̅-l̅0)-34(sin⁡2l̅-sin⁡2l̅0)]λk2},
with
(2.23)α(t)=α0+4a̅5μ3{(l̅-l̅0)λα+(sin⁡l̅-sin⁡l̅0)λh-(cos⁡l̅-cos⁡l̅0)λk},(2.24)h(t)=h0+a̅5μ3{4(sin⁡l̅-sin⁡l̅0)λα+[52(l̅-l̅0)+34(sin⁡2l̅-sin⁡2l̅0)]λh-34(cos⁡2l̅-cos⁡2l̅0)λk},(2.25)k(t)=k0+a̅5μ3{-4(cos⁡l̅-cos⁡l̅0)λα-34(cos⁡2l̅-cos⁡2l̅0)λh+[52(l̅-l̅0)-34(sin⁡2l̅-sin⁡2l̅0)]λk},(2.26)λα=12μ3a̅5{Δα(5Δl̅+3sin⁡Δl̅)10Δl̅2+6Δl̅sin⁡Δl̅-64sin2(Δl̅/2)},(2.27)λh=-μ3a̅5{8Δαsin⁡Δl̅/210Δl̅2+6Δl̅sin⁡Δl̅-64sin2(Δl̅/2)},(2.28)λk=0,
where α=a/a̅,h=ecosω,k=esinω, where a is the semimajor axis, e is the eccentricity, ω is the argument of the pericenter, ℓ̅=ℓ̅0+n̅(t-t0), 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:
(2.29)Γ*=1n̅a̅{(λhsin⁡l̅-λkcos⁡l̅)er+2(λα+λhcos⁡l̅+λksin⁡l̅)es},
where er and es 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̅=(a0+af)/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.


## 2.1. 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/g0, is between 10−4 and 10−2. For such system, the fuel consumption is described by the variable J defined as
(2.1)J=12∫t0tγ2dt,
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:
(2.2)J=Pmax⁡(1m-1m0),
where Pmax is the maximum power, and m0 is the initial mass. The minimization of the final value Jf is equivalent to the maximization of mf or the minimization of the fuel consumption.The optimization problem concerning with simple transfers (no rendezvous) between coplanar orbits is formulated as: at timet, 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]:(2.3)dudt=v2r-μr2+R,dvdt=-uvr+S,drdt=u,dJdt=12(R2+S2),
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 t0=0:
(2.4)u(0)=0v(0)=1r(0)=1J(0)=0,
to the final state at the prescribed final time tf:
(2.5)u(tf)=0v(tf)=μrfr(tf)=rf,
such that Jf is a minimum, that is, the performance index is defined by:(2.6)IP=J(tf).
Equations (2.4) and (2.5) are given in canonical units, and they define the initial and final circular orbits. u(tf),v(tf), and r(tf) denote the state variables at the prescribed final time tf, and 0, μ/rf, and rf are the prescribed values defining the final circular orbit. Similar definition applies at the initial time t0=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.


## 2.2. 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):(2.7)IP=J(tf)+k1(u(tf))2+k2(v(tf)-1rf)2+k3(r(tf)-rf)2,
where k1,k2,k3 ≫ 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]:(2.8)H=λu(v2r-μr2+R)+λv(-uvr+S)+λru+12λJ(R2+S2).From the HamiltonianH, one finds the adjoint equations:
(2.9)dλudt=vrλv-λr,dλvdt=-2vrλu+urλv,dλrdt=(v2r2-2μr3)λu-uvr2λv,dλJdt=0,
and, from the performance index defined by (2.7), we get the terminal conditions for the adjoint equations:
(2.10)λu(tf)=-2k1u(tf),λv(tf)=-2k2(v(tf)-1rf),λr(tf)=-2k3(r(tf)-rf),λJ(tf)=-1.The algorithm also requires the partial derivatives of the HamiltonianH with respect to the control variables. These partial derivatives are given by:(2.11)∂H∂R=λu+RλJ,∂H∂S=λv+SλJ.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:
(2.12)Hαα=[λJ00λJ],Hλα=[100100RS],Hxα=0(4×2nullmatrix),Hλx=[02vr-v2r2+2μr30-vr-uruvr2010000000],Hxx=[0-λvrvr2λv0-λvr2λur-2vr2λu+ur2λv0λvvr2-2vr2λu+ur2λv(2v2r3-6μr4)λu-2uvr3λv00000].α denotes the control vector αT=[RS], and it has been introduced to avoid confusion with the state variable u,x is the state vector xT=[uvrJ], and λ is the adjoint vector λT=[λuλvλrλJ].From (2.5) and (2.6), one defines the functions ψ  and Φ:(2.13)ψ=[u(tf)v(tf)-1rfr(tf)-rf],(2.14)Φ=J(tf)+Λ1u(tf)+Λ2(v(tf)-1rf)+Λ3(r(tf)-rf),
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)ψx=[100001000010],Φxx=0.(4×4nullmatrix).The results of the gradient-based algorithm to the optimization problem described above are presented in Section3.


## 2.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,
(2.16)R*=-λuλJ,S*=-λvλJ.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):
(2.17)dudt=v2r-μr2-λuλJ,dvdt=-uvr-λvλJ,drdt=u,dJdt=12λJ2(λu2+λv2),dλudt=vrλv-λr,dλvdt=-2vrλu+urλv,dλrdt=(v2r2-2μr3)λu-uvr2λv,dλJdt=0,
with the boundary conditions:
(2.18)u(0)=0,v(0)=1,r(0)=1,J(0)=0,u(tf)=0,v(tf)=μrf,r(tf)=rf,λJ(tf)=-1.The neighboring extremals algorithms are based on the solution of a linearized two-point boundary value problem that involves the derivatives of the right-hand side of (2.17) with respect to the state and adjoint variables [16, 21, 27]. These equations can be put in the following form:
(2.19)dδxdt=Aδx+Bδλ,dδλdt=Cδx-ATδλ,
with δx(t)=xn+1(t)-xn(t) and δλ(t)=λn+1(t)-λn(t), where n denotes the iterate, and A, B, and C are matrices given by:
(2.20)A=[02vr-a0-vr-uruvr2010000000],B=[-1λJ00λuλJ20-1λJ0λvλJ20000λuλJ2λvλJ20-cλJ3],C=[0λvr-vλvr20λvr-2λurb0-vλvr2bd00000],
where
(2.21)a=v2r2-2μr3,b=2vr2λu-ur2λv,c=λu2+λv2,d=(-2v2r3+6μr4)λu+2uvr3λv.The results of the neighboring extremals algorithms to the optimization problem described above are presented in Section3.


## 2.4. 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:
(2.22)u=na(hsin⁡l̅-kcos⁡l̅),v=na(1+hcos⁡l̅+ksin⁡l̅),r=a1+hcos⁡l̅+ksin⁡l̅,J(t)=12a̅5μ3{[52(l̅-l̅0)-34(sin⁡2l̅-sin⁡2l̅0)]λk24(l̅-l̅0)λα2+8(sin⁡l̅-sin⁡l̅0)λαλh-8(cos⁡l̅-cos⁡l̅0)λαλk+[52(l̅-l̅0)+34(sin⁡2l̅-sin⁡2l̅0)]λh2-32(cos⁡2l̅-cos⁡2l̅0)λhλk+[52(l̅-l̅0)-34(sin⁡2l̅-sin⁡2l̅0)]λk2},
with
(2.23)α(t)=α0+4a̅5μ3{(l̅-l̅0)λα+(sin⁡l̅-sin⁡l̅0)λh-(cos⁡l̅-cos⁡l̅0)λk},(2.24)h(t)=h0+a̅5μ3{4(sin⁡l̅-sin⁡l̅0)λα+[52(l̅-l̅0)+34(sin⁡2l̅-sin⁡2l̅0)]λh-34(cos⁡2l̅-cos⁡2l̅0)λk},(2.25)k(t)=k0+a̅5μ3{-4(cos⁡l̅-cos⁡l̅0)λα-34(cos⁡2l̅-cos⁡2l̅0)λh+[52(l̅-l̅0)-34(sin⁡2l̅-sin⁡2l̅0)]λk},(2.26)λα=12μ3a̅5{Δα(5Δl̅+3sin⁡Δl̅)10Δl̅2+6Δl̅sin⁡Δl̅-64sin2(Δl̅/2)},(2.27)λh=-μ3a̅5{8Δαsin⁡Δl̅/210Δl̅2+6Δl̅sin⁡Δl̅-64sin2(Δl̅/2)},(2.28)λk=0,
where α=a/a̅,h=ecosω,k=esinω, where a is the semimajor axis, e is the eccentricity, ω is the argument of the pericenter, ℓ̅=ℓ̅0+n̅(t-t0), 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:
(2.29)Γ*=1n̅a̅{(λhsin⁡l̅-λkcos⁡l̅)er+2(λα+λhcos⁡l̅+λksin⁡l̅)es},
where er and es 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̅=(a0+af)/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.


## 3. 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ρ=rf/r0 and for various time of flight tf-t0 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.Table 1
Consumption variableJ(ρ>1) for transfers with small time of flight.
ρtf-t0JanalJgradJNeigh1JNeigh2drel1drel2drel31.02502.0 3.0 4.0 5.03.5856 × 10‒48.4459 × 10‒53.1226 × 10‒51.7138 × 10‒53.5855 × 10‒48.4462 × 10‒53.1233 × 10‒51.7147 × 10‒53.5854 × 10‒48.4456 × 10‒53.1230 × 10‒51.7143 × 10‒53.5854 × 10‒48.4456 × 10‒53.1230 × 10‒51.7143 × 10‒50.000.00 0.01 0.030.000.01 0.010.020.000.00 0.00 0.001.05002.03.0 4.0 5.01.4463 × 10‒33.4169 × 10‒41.2533 × 10‒46.7541 × 10‒51.4463 × 10‒33.4166 × 10‒41.2538 × 10‒46.7611 × 10‒51.4459 × 10‒33.4164 × 10−41.2537 × 10−46.7598 × 10−51.4459 × 10−33.4164 × 10−41.2537 × 10−46.7598 × 10−50.030.01 0.03 0.080.030.01 0.00 0.020.000.00 0.00 0.001.10002.03.0 4.0 5.05.8778 × 10−31.3977 × 10−35.0619 × 10−42.6374 × 10−45.8741 × 10−31.3970 × 10−35.0666 × 10−42.6453 × 10−45.8716 × 10−31.3969 × 10−35.0664 × 10−42.6451 × 10−45.8716 × 10−31.3969 × 10−35.0664 × 10−42.6451 × 10−40.110.06 0.09  0.290.040.00 0.00 0.010.000.00 0.00 0.001.20002.03.0 4.0 5.02.4187 × 10−25.8370 × 10−32.0813 × 10−31.0260 × 10−32.4097 × 10−25.8200 × 10−32.0845 × 10−31.0346 × 10−32.4097 × 10−25.8199 × 10−32.0844 × 10−31.0345 × 10−32.4097 × 10−25.8199 × 10−32.0844 × 10−31.0345 × 10−30.370.29 0.15 0.820.000.00 0.00 0.000.000.00 0.00  0.00Table 2
Consumption variableJ(ρ<1) for transfers with small time of flight.
ρtf-t0JlinearJgradJNeigh1JNeigh2drel1drel2drel30.80002.03.0 4.0 5.02.0951 × 10−24.9040 × 10−32.0703 × 10−31.3838 × 10−32.0842 × 10−24.9173 × 10−32.1047 × 10−31.4198 × 10−32.0842 × 10−24.9172 × 10−32.1046 × 10−31.4197 × 10−32.0842 × 10−24.9172 × 10−32.1046 × 10−31.4197 × 10−30.520.27 1.63 2.530.000.00 0.00 0.000.000.00 0.00 0.000.90002.03.0 4.0 5.05.4740 × 10−31.2771 × 10−35.0063 × 10−43.0496 × 10−45.4672 × 10−31.2772 × 10−35.0198 × 10−43.0653 × 10−45.4671 × 10−31.2771 × 10−35.0198 × 10−43.0652 × 10−45.4671 × 10−31.2771 × 10−35.0198 × 10−43.0652 × 10−40.130.00 0.27 0.510.000.01 0.00 0.000.000.00 0.00 0.000.95002.03.04.0 5.01.3958 × 10‒33.2649 × 10−41.2451 × 10−47.2585 × 10−51.3955 × 10−33.2649 × 10−41.2459 × 10−47.2671 × 10−51.3955 × 10−33.2647 × 10−41.2458 × 10−47.2667 × 10−51.3955 × 10−33.2647 × 10−41.2458 × 10−47.2667 × 10−50.020.01 0.06 0.110.000.01 0.01 0.010.000.00 0.00 0.000.97502.03.0 4.0 5.03.5225 × 10−48.2555 × 10−53.1120 × 10−51.7765 × 10−53.5231 × 10−48.2560 × 10−53.1126 × 10−51.7772 × 10−53.5223 × 10−48.2554 × 10−53.1124 × 10−51.7771 × 10−53.5223 × 10−48.2553 × 10−53.1124 × 10−51.7771 × 10−50.010.00 0.01 0.030.020.01 0.01 0.000.000.00 0.00 0.00Table 3
Consumption variableJ(ρ>1) for transfers with moderate time of flight.
ρtf-t0JlinearJgradJNeigh1JNeigh2drel1drel2drel31.025020.030.0 40.0 50.03.7722 × 10−62.5221 × 10−61.8855 × 10−61.5066 × 10−63.7733 × 10−62.5226 × 10−61.8859 × 10−61.5072 × 10−63.7730 × 10−62.5230 × 10−61.8860 × 10−61.5070 × 10−63.7714 × 10−62.5209 × 10−61.8841 × 10−61.5049 × 10−60.020.04 0.03 0.030.010.02 0.00 0.010.040.08 0.10 0.141.050020.030.0 40.0 50.01.4520 × 10−59.7411 × 10−67.2599 × 10−65.8158 × 10−61.4536 × 10−59.7482 × 10−67.2659 × 10−65.8199 × 10−61.4533 × 10−59.7480 × 10−67.2660 × 10−65.8200 × 10−61.4531 × 10−59.7460 × 10−67.2636 × 10−65.8182 × 10−60.090.07 0.08 0.070.020.00 0.00 0.000.010.02 0.03 0.031.100020.030.0 40.0 50.05.4007 × 10−53.6278 × 10−52.7003 × 10−52.1653 × 10−55.4168 × 10−53.6390 × 10−52.7083 × 10−52.1719 × 10−55.4167 × 10−53.6389 × 10−52.7078 × 10−52.1718 × 10−55.4165 × 10−53.6387 × 10−52.7077 × 10−52.1716 × 10−50.300.30 0.27 0.300.000.00 0.02 0.000.000.01 0.00 0.011.200020.030.0 40.0 50.01.8980 × 10−41.2543 × 10−49.4416 × 10−57.5157 × 10−51.9172 × 10−41.2695 × 10−49.5396 × 10−57.5976 × 10−51.9154 × 10−41.2693 × 10−49.5391 × 10−57.5928 × 10−51.9154 × 10−41.2693 × 10−49.5390 × 10−57.5927 × 10−50.911.18 1.02 1.020.090.02 0.01 0.060.000.00 0.00 0.00Table 4
Consumption variableJ(ρ<1) for transfers with moderate time of flight.
ρtf-t0JlinearJgradJNeigh1JNeigh2drel1drel2drel30.800020.0 30.0 40.0 50.03.4529 × 10−4 2.2973 × 10−41.7196 × 10−4 1.3736 × 10−43.5015 × 10−4 2.3313 × 10−41.7467 × 10−41.3978 × 10−43.5011 × 10−42.3316 × 10−41.7465 × 10−4 1.3959 × 10−43.5010 × 10−42.3311 × 10−4 1.7465 × 10−41.3959 × 10−41.37 1.47 1.54 1.590.010.01 0.01 0.130.000.02 0.00 0.000.900020.030.0 40.0 50.07.3862 × 10−54.8663 × 10−53.6467 × 10−5 2.9218 × 10−57.4146 × 10−54.8851 × 10−53.6588 × 10−5 2.9317 × 10−57.4146 × 10−54.8852 × 10−53.6589 × 10−52.9316 × 10−57.4144 × 10−54.8850 × 10−53.6587 × 10−52.9314 × 10−50.380.39 0.33 0.330.000.00 0.00 0.000.000.00 0.00 0.010.950020.030.0 40.0 50.01.7023 × 10−5 1.1240 × 10−58.4569 × 10−6 6.7519 × 10−61.7042 × 10−5 1.1251 × 10−58.4642 × 10−6 6.7581 × 10−61.7040 × 10−51.1249 × 10−58.4640 × 10−66.7580 × 10−61.7038 × 10−51.1247 × 10−5 8.4620 × 10−66.7561 × 10−60.100.08 0.08 0.090.010.02 0.00 0.000.010.02 0.02 0.030.975020.030.0 40.050.04.0858 × 10−6 2.7075 × 10−62.0361 × 10−6 1.6230 × 10−64.0869 × 10−62.7081 × 10−6 2.0366 × 10−6 1.6234 × 10−64.0870 × 10−62.7080 × 10−62.0360 × 10−6 1.6230 × 10−64.0847 × 10−6 2.7059 × 10−62.0349 × 10−6 1.6216 × 10−60.03 0.02 0.00 0.000.000.00 0.03 0.020.060.08 0.05 0.09Table 5
Consumption variableJ for Earth-Venus transfers.
ρtf-t0JlinearJgradJNeigh1JNeigh2drel1drel2drel30.72702.03.0 4.0 5.03.7654 × 10−28.9269 × 10−34.0482 × 10−32.8941 × 10−33.7299 × 10−29.0261 × 10−3 4.2133 × 10−33.0573 × 10−33.7298 × 10−2 9.0259 × 10−34.2131 × 10−33.0572 × 10−33.7298 × 10−29.0259 × 10−34.2131 × 10−33.0572 × 10−30.951.10 3.915.330.00 0.00 0.00 0.000.00 0.00 0.00 0.0020.030.0 40.0 50.07.2355 × 10−4 4.8236 × 10−43.6177 × 10−42.8942 × 10−47.4857 × 10−44.9863 × 10−43.7385 × 10−42.9903 × 10−47.4856 × 10−44.9862 × 10−43.7384 × 10−42.9901 × 10−47.4856 × 10−44.9862 × 10−43.7383 × 10−4 2.9897 × 10−43.343.26 3.22 3.200.00 0.00 0.00 0.010.000.00 0.00 0.01Table 6
Consumption variableJ for Earth-Mars transfers.
ρtf-t0JlinearJgradJNeigh1JNeigh2drel1drel2drel31.52362.03.0 4.0 5.01.7743 × 10−14.4947 × 10−21.6051 × 10−2 7.2498 × 10−31.7434 × 10−14.4067 × 10−21.5889 × 10−27.3352 × 10−31.7434 × 10−1 4.4066 × 10−2 1.5889 × 10−27.3351 × 10−31.7434 × 10−1 4.4066 × 10−21.5889 × 10−27.3351 × 10−31.771.99 1.02 1.160.000.00 0.00 0.000.00 0.00 0.00 0.0020.030.0 40.0 50.08.6591 × 10−4 5.7537 × 10−44.2991 × 10−4 3.4273 × 10−49.3232 × 10−4 6.1074 × 10−44.5311 × 10−4 3.6096 × 10−49.3151 × 10−46.1071 × 10−44.5296 × 10−4 3.6093 × 10−49.3158 × 10−4 6.1073 × 10−44.5299 × 10−43.6095 × 10−47.025.79 5.09 5.040.09 0.00 0.03 0.010.010.00 0.01 0.01Table 7
Consumption variableJ for interplanetary transfers with large-amplitude.
ρtf-t0JNeigh1JNeigh22.50020.030.0 40.0 50.0 60.03.7736 × 10−32.3900 × 10−31.7541 × 10−31.3859 × 10−31.1454 × 10−33.7733 × 10−32.3900 × 10−31.7541 × 10−31.3859 × 10−31.1453 × 10−35.20320.030.0 40.0 50.0 60.01.3746 × 10−27.5307 × 10−35.0533 × 10−33.7103 × 10−32.9897 × 10−31.3746 × 10−27.5309 × 10−3 5.0533 × 10−33.7100 × 10−32.9896 × 10−39.51960.070.0 80.0 90.0 100.05.9485 × 10−34.6970 × 10−33.9003 × 10−33.3295 × 10−32.8858 × 10−35.9392 × 10−34.6980 × 10−33.9009 × 10−33.3285 × 10−32.8857 × 10−3Table 8
Consumption variableJ for large transfers.
ρtf-t0JNeigh1JNeigh22.500100.0 125.0 150.0 175.0 200.06.7827 × 10−45.4241 × 10−44.5154 × 10−43.8651 × 10−43.3812 × 10−46.7818 × 10−4 5.4241 × 10−44.5154 × 10−43.8688 × 10−43.3812 × 10−43.750100.0125.0 150.0 175.0 200.01.1975 × 10−39.5018 × 10−47.8787 × 10−46.7324 × 10−45.8793 × 10−41.1966 × 10−39.4987 × 10−47.8784 × 10−46.7325 × 10−45.8794 × 10−45.000100.0125.0 150.0 175.0 200.01.6105 × 10−31.2678 × 10−31.0441 × 10−38.8948 × 10−47.7584 × 10−41.6106 × 10−31.2678 × 10−31.0441 × 10−38.8947 × 10−47.7584 × 10−4Tables1–4 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:(3.1)drel⁡1=|(Jneigh1-Jlinear)Jneigh1|×100%,drel⁡2=|(Jneigh1-Jgrad)Jneigh1|×100%,drel⁡3=|(Jneigh1-Jneigh2)Jneigh1|×100%.
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 Tables5, 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.From the results presented in Tables1–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, drel1<0.5%;(2)
For transfers with small time of flight (tf-t0=2.0,3.0,4.0,5.0 time units), Tables 1, 2, 5 and 6 show that the maximum absolute relative difference drel1 occur for the most of the transfers with tf-t0=5. This maximum value of drel1 is about 2% for ρ>1 and 5.5% for ρ<1;(3)
For transfers with moderate time of flight (tf-t0=20.0,30.0,40.0,50.0 time units), Tables 3, 4, 5, and 6 show that the maximum absolute relative difference drel1 is about 7% for ρ>1, and 3.5% for ρ<1;(4)
In all cases described above, the maximum absolute relative differencesdrel1 occur for transfers with large radial excursion;(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 flighttf-t0=2.0 (time units) is approximately ten times the fuel consumption for transfers with time of flight tf-t0=4.0 (time units), and, it is approximately hundred times the fuel consumption for transfers with time of flight tf-t0=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%>drel1>1.0%;(8)
Tables1–6 show that the results obtained through the numerical algorithms—gradient and neighboring extremals—are very quite similar, regardless the amplitude of maneuver and the time of flight;(9)
Table7 shows that an Earth-asteroid belt mission with tf-t0=20.0 time units (approximately, 3.2 years) and an Earth-Jupiter mission with tf-t0=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 tf-t0 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 tf-t0=30.0 time units (approximately, 4.8 years) and an Earth-Venus mission with tf-t0=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. Figures1, 2, 3, and 4 illustrate these plots for ρ= 0.727, 0.950, 0.975, 1.025, 1.050, and 1.523, with tf-t0=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.Figure 1
Thrust acceleration fortf-t0=2.0.Figure 2
Thrust acceleration fortf-t0=3.0.Figure 3
Thrust acceleration fortf-t0=30.0.Figure 4
Thrust acceleration fortf-t0=50.0.Similarly, Figures5 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.Figure 5
Thrust acceleration for large-amplitude transfers withtf-t0=100.0.Figure 6
Thrust acceleration for large-amplitude transfers withtf-t0=200.0.Figure 7
Consumption for Earth-Venus and Earth-Mars transfers with moderate time of flight.Figure 8
Consumption for Earth-asteroid belt and Earth-Jupiter transfers.Figure 9
Consumption for large transfers.Figures3, 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.Figures7 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 tf-t0 can be performed with the same amount of fuel (see Comment  9).


## 4. 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.



---

*Source: 168632-2012-05-27.xml*