# A First-Order Analytical Theory for Optimal Low-Thrust Limited-Power Transfers between Arbitrary Elliptical Coplanar Orbits

Authors: Sandro da Silva Fernandes; Francisco das Chagas Carvalho
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

A complete first-order analytical solution, which includes the short periodic terms, for the problem of optimal low-thrust limited-power transfers between arbitrary elliptic coplanar orbits in a Newtonian central gravity field is obtained through canonical transformation theory. The optimization problem is formulated as a Mayer problem of optimal control theory with Cartesian elements—position and velocity vectors—as state variables. After applying the Pontryagin maximum principle and determining the maximum Hamiltonian, classical orbital elements are introduced through a Mathieu transformation. The short periodic terms are then eliminated from the maximum Hamiltonian through an infinitesimal canonical transformation built through Hori method. Closed-form analytical solutions are obtained for the average canonical system by solving the Hamilton-Jacobi equation through separation of variables technique. For transfers between close orbits a simplified solution is straightforwardly derived by linearizing the new Hamiltonian and the generating function obtained through Hori method.


---

## Body

## 1. Introduction

In the general analysis of optimal
space trajectories, two idealized propulsion models are of most frequent use
[1]: LP and CEV systems. In the power-limited variable ejection velocity
systems—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—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 about more realistic
problems.The purpose of this work is to
present a complete first-order analytical solution, which includes the short periodic
terms, for the problem of optimal low-thrust limited-power transfers (no
rendezvous) between arbitrary elliptic coplanar orbits in a Newtonian central
gravity field. This 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.The
optimization problem is formulated as a Mayer problem of optimal control theory
with Cartesian elements—position and
velocity vectors—as state
variables [1]. After applying the Pontryagin Maximum Principle and determining
the maximum Hamiltonian, classical orbital elements are introduced through a
canonical Mathieu transformation. The short periodic terms are then eliminated
from the maximum Hamiltonian through an infinitesimal canonical transformation
built through Hori method—a perturbation
canonical method based on Lie series [2]. The new maximum Hamiltonian,
resulting from the infinitesimal canonical transformation, describes the
extremal trajectories associated with the long duration maneuvers.Closed-form
analytical solutions are obtained for the average canonical system governed by
this new Hamiltonian. The separation of variables technique is applied to solve
the Hamilton-Jacobi equation associated to the average canonical system [3]. For long duration maneuvers, the existence of conjugate points is
investigated through the Jacobi condition and the envelope of extremals is
obtained considering different configurations of initial and final orbits.
Curves representing maneuvers with same fuel consumption for a specified
duration, called isocost curves, are presented for several configurations of
initial and final orbits. 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.For
transfers between close coplanar elliptic orbits, a simplified solution is
straightforwardly derived by linearizing the new Hamiltonian and the generating
function obtained through Hori method. Since this solution is given by a linear
system of algebraic equations in the initial value of the adjoint variables,
the two-point boundary value problem is solved by simple techniques.Some
of the results presented in the paper are very similar to the ones obtained by
Edelbaum [4] and Marec and Vinh [5] through different approaches based on the
concept of mean Hamiltonian. We note that the theory described in the paper
includes directly the periodic terms such that it can be applied to any
transfer independently of its duration. On the other
hand, since the theory is based on canonical transformations, a first-order
solution for the Cartesian elements can also be easily obtained from Lie
theorem [2].


## 2. Optimal Space Trajectories

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/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∫t0tfγ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 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 of the fuel consumption Jf is equivalent to the maximization of mf.Let us consider the motion of a space vehicleM, powered by a limited-power engine in
a Newtonian central gravity field. At time t,
the state of the vehicle is defined by the position vector r(t),
the velocity vector v(t) and the consumption variable J(t).
This set of variables provides a straightforward physical interpretation to the
optimal thrust acceleration.The problem of general
optimal low-thrust limited-power transfer can be formulated as a Mayer problem
of optimal control as follows [1]. It is proposed to transfer the space vehicle M from the initial state (r0,v0,0) at time t0 to the final state (rf,vf,Jf) at time tf,
such that the final consumption variable Jf is a minimum (Figure 1).Figure 1
Geometry of the transfer problem.For coplanar transfer
problem, the initial and final conditions are specified in terms of orbital
elements—semimajor axis,
eccentricity, and pericenter argument—introduced in the
next section. It is also assumed thattf−t0 and the position of the vehicle in the initial
orbit are specified. The state equations are (2.3)drdt=v,dvdt=−μr3r+γ,dJdt=12γ2,
where μ is the gravitational parameter. It is assumed
that the thrust direction and magnitude are unconstrained [1].According to the
Pontryagin Maximum Principle [6], the optimal thrust acceleration γ* must be selected from the admissible controls
such that the Hamiltonian H reaches
its maximum. The Hamiltonian is formed using (2.3), (2.4)H=pr•v+pv•(−μr3r+γ)+12pJγ2,where pr, pv, and pJ are the adjoint variables and dot denotes the dot
product. Since the optimization problem is unconstrained, it follows that (2.5)γ*=−pvpJ. The optimal thrust
acceleration γ* is modulated [1]. Then, the optimal
trajectories are governed by the maximum Hamiltonian H*, (2.6)H*=pr•v−pv•μr3r−pv22pJ.Since J is an ignorable variable, it is
determined by simple integration, and its adjoint pJ is a first integral, whose value is obtained
from the transversality conditions: pJ=−1.
Thus, (2.7)H*=pr•v−pv•μr3r+pv22.This Hamiltonian can be decomposed into two parts:(2.8)H0=pr•v−pv•μr3r,Hγ*=pv22,where H0 is the undisturbed Hamiltonian and Hγ* is the disturbing function concerning the
optimal thrust acceleration. The undisturbed Hamiltonian H0 plays a fundamental role in the canonical
transformation theory described in the next section.


## 3. Transformation from Cartesian Elements to Orbital Elements

In
order to obtain a first-order analytical solution for transfers between arbitrary
elliptic coplanar orbits, the set of classical orbital elements is introduced
through a canonical transformation by using the properties of generalized
canonical systems [7].Consider
the extended two-body problem defined by the canonical system of differential
equations governed by the undisturbed HamiltonianH0,
that is, the problem defined by the differential equations for the state
variables r and v,
which describe the classical two-body problem, and the differential equations
for the adjoint variables pr and pv.
These equations are given by (3.1)drdt=vdprdt=μr3(pv−3(pv•er)er),dvdt=−μr3rdpvdt=−pr,    where er is the unit vector pointing radially outward
of the moving frame of reference.The general solution of the differential equations for
the state variablesr and v is well-known from the classical
two-body problem [8]. In the coplanar case, it is given by (3.2)r=a(1−e2)1+ecosfer,(3.3)v=μa(1−e2)[esinfer+(1+ecosf)es],where es is the
unit vector along circumferential direction of the moving frame of reference; a is the semimajor axis, e is the eccentricity, and f is the true anomaly. The unit vectors er and es are
written as function of the orbital elements as follows: (3.4)er=cos(ω+f)i+sin(ω+f)j,es=−sin(ω+f)i+cos(ω+f)j where i and j are the unit vectors of
a fixed frame of reference (Figure 1) and ω is the pericenter argument.The
general solution of the differential equations for the adjoint variablespr and pv is obtained
by computing the inverse of the Jacobian matrix of the point transformation
between the Cartesian elements and the orbital ones, defined by (3.2) and (3.3). This
matrix is obtained through the variations of the orbital elements induced by the
variations in the Cartesian elements, as described in the next paragraphs.Let
us consider the inverse of the point transformation defined by (3.2) and (3.3): (3.5)a=r2−(rv2/μ),(3.6)e2=1−h2μa,(3.7)cosω=i•ee,(3.8)cosE=1e(1−ra),where the eccentricity vector e and the angular momentum vector h, shown in Figure 2, are given, as
function of the Cartesian elements r and v, by the following equations: (3.9)e=1μ[(v2−μr)r−(r•v)v],h=r×v.
Here, the symbol × denotes the cross product. The true anomaly f has been replaced by eccentric anomaly E. These anomalies are related
through the equation (3.10)tanf2=1+e1−etanE2.Figure 2
The vectorse and h.Now,
consider the variations in the Cartesian elements,δr and δv,
given in the moving frame of reference by (3.11)δr=δξer+δηes,δv=δuer+δves.The
variations of the orbital elements—a, e, ω, and E—induced by the
variations in the Cartesian elements, δr and δv,
are obtained straightforwardly from (3.5) through (3.8) and are given by (3.12)δa=2(ar)2r•δrr+2a2μv•δv,δe=−h•δhμae+h22μa2eδa,δω=δeecotω−i•δeesinω,δE=1esinE[r•δrra−ra2δa+cosEδe],
where the variations of the vectors e and h are written as (3.13)μδe=(vs2δξ−vsvrδη+2rvsδv)er+(−vsvrδξ+(vr2−μr)δη−rvsδu−rvrδv)es,(3.14)δh=(vsδξ−vrδη+rδv)ew,being the unit vector ew normal
to the plane of the orbit. Here, vr and vs denote
the radial and circumferential components of the velocity vector, respectively
(see (3.3)).From
(3.2), (3.3), and (3.9) through (3.14), one gets the explicit form of the variations
of the orbital elements—a, e, ω, and E—induced by the
variations in the Cartesian elements, δr and δv.
At this point, it is useful to replace the variation in the eccentric anomaly E by the variation in the mean anomaly M, obtained from the well-known Kepler
equation (3.15)M=E−esinE.This variation is given by
(3.16)δM=(ra)δE−sinEδe.Thus, (3.17)δa=2(ar)2δξ+2esinfn1−e2δu+21−e2n(ar)δv,δe=a(1−e2)r2cosEδξ+sinfaδη+1−e2nasinfδu+1−e2na(cosE+cosf)δv,δω=sinferδξ−e+cosfae(1−e2)δη−1−e2naecosfδu+1−e2naesinf(1+11+ecosf)δv,δM=−1−e3cosEer1−e2sinfδξ+1−e2aecosfδη+(1−e2)nae(cosf−2e1+ecosf)δu−(1−e2)naesinf(1+11+ecosf)δv,
where  n=μ/a3 is the
mean motion.Equation
(3.17) can be put in the form (3.18)[δaδeδωδM]=Δ−1[δξδηδuδv],where the matrix Δ−1 is inverse
Jacobian matrix of the point transformation between the Cartesian elements and
the orbital ones, defined by (3.2) and (3.3).Following
the properties of generalized canonical systems [7], the general solution of
the differential equations for the adjoint variables pr and pv is given
by (3.19)[prpv]=(Δ−1)T[papepωpM],with pr and pv expressed in the moving frame of reference.
Thus, it follows from (3.17) through (3.19) that
(3.20)pr=ar2{2apa+((1−e2)cosE)pe+(ra)sinfe(pω−(1−e3cosE)1−e2pM)}er+{sinfape−(e+cosf)ae(1−e2)pω+1−e2cosfaepM}es,pv=1na1−e2{{2aesinfpa+((1−e2)sinf)pe−(1−e2)cosfepω+(1−e2)3/2e(cosf−2e1+ecosf)pM}er+{2a(1−e2)(ar)pa+(1−e2)(cosf+cosE)pe+(1−e2)sinfe(1+11+ecosf)(pω−1−e2pM)}es}.(r/a), (r/a)sinf, and so forth are known functions of the
elliptic motion, which can be written explicitly in terms of eccentricity and
mean anomaly through Lagrange series [8].Equations (3.2), (3.3), and (3.20) define a
canonical Mathieu transformation between the Cartesian elements (r,v,pr,pv) and the
orbital ones (a,e,ω,M,pa,pe,pω,pM). The Hamiltonian H* is
invariant with respect to this canonical transformation. Thus, (3.21)H0=npM,(3.22)Hγ*=12n2a2(1−e2)×{12(1−cos2f)[2aepa+(1−e2)pe]2+2(1−e2)sin2f[−apapω−(1−e2)2epepω]+4(1−e2)3/2sinf(−2e1+ecosf+cosf)[apapM+(1−e2)2epepM]+(1−e2)22e2(1+cos2f)pω2−2(1−e2)5/2e2(−2e1+ecosf+cosf)cosfpωpM+(1−e2)3e2(−2e1+ecosf+cosf)2pM2+4a2(1−e2)2(ar)2pa2+4a(1−e2)2(ar)(cosE+cosf)pape+(1−e2)2(cosE+cosf)2pe2+4a(1−e2)2e(ar)sinf(1+11+ecosf)[papω−1−e2papM]+2(1−e2)2e(cosE+cosf)(1+11+ecosf)sinf[pepω−1−e2pepM]+[(1−e2)e(1+11+ecosf)sinf[pω−1−e2pM]]2}.The new Hamiltonian functionH*,
defined by (3.21) and (3.22), describes the optimal low-thrust limited-power
trajectories in a Newtonian central gravity field for transfers between
arbitrary elliptic coplanar orbits. It should be noted that H* becomes singular for circular orbits and,
thus, the theory developed in the paper cannot be applied to transfers
involving orbits with very small eccentricities.A first-order analytical solution for the canonical
system described by the Hamiltonian functionH* will be derived in the next sections through
canonical transformation theory.


## 4. A First-Order Analytical Solution

In order to get an approximate formal solution to
the problem of transfers between arbitrary elliptic coplanar orbits, a
perturbation technique based on Lie series—Hori method—will be applied. For
completeness, a brief description of Hori method [2] is presented in the
appendix.In
this study, for low-thrust propulsion systems, it is assumed that the functionsH0 and Hγ* have
different order of magnitude: H0 has
zero order and Hγ* has
first order in a small parameter associated to the magnitude of the thrust
acceleration.Consider
an infinitesimal canonical transformation built through Hori method,(4.1)(a,e,ω,M,pa,pe,pω,pM)→(a′,e′,ω′,M′,pa′,pe′,pω′,pM′). Following the algorithm described in the
appendix, at order 0 one gets (4.2)F0=n′pM′.Now,
consider the undisturbed canonical system described byF0 : (4.3)da′dt=0,dpa′dt=32n′a′pM′,de′dt=0,dpe′dt=0,dω′dt=0,dpω′dt=0,dM′dt=n,dpM′dt=0,    general solution of which is given by (4.4)a′=a0′,e′=e0′,ω′=ω0′,M′=M0′+n′(t−t0),pa′=pa0′+32n′(t−t0)a′pM′,pe′=pe0′,pω′=pω0′,pM′=pM0′.The subscript 0 denotes the constants of
integration.Introducing
the general solution defined by (4.4) in the equation of order 1 of the algorithm
(A.6), one gets (4.5)∂S1∂t=Hγ*−F1.The functions Hγ* and F1 are
written in terms of the constants of integration. According to (A.11), the mean value of Hγ* must be calculated.
After lengthy calculations, one gets (4.6)F1=a′2μ{4a′2p′a2+52(1−e′2)p′e2+(5−4e′2)2e′2p′ω2},(4.7)S1=12a′5μ3{8e′sinE′a′2pa′2+8(1−e′2)sinE′a′pa′pe′−8(1−e′2)1/2e′cosE′a′pa′pω′+(1−e′2)[−54e′sinE′+34sin2E′−112e′sin3E′]pe′2+2(1−e′2)1/2e′[54e′cosE′+14(e′2−3)cos2E′+112e′cos3E′]pe′pω′+1e′2[(54−e′2)e′sinE′−12(32−e′2)sin2E′+112e′sin3E′]pω′2}.
Terms factored by pM′ have
been omitted in equations above, since only transfers (no rendezvous) are
considered [4]. But, it should be noted that mixed secular terms can arise in
the new Hamiltonian in the case of rendezvous maneuvers. For a first-order
solution, F′=F1 and S=S1.The new HamiltonianF′ governs the optimal long duration transfers
between arbitrary elliptic coplanar orbits. The general solution of the
canonical system described by this Hamiltonian will be obtained through two
canonical transformations.First, consider the Mathieu transformation:(4.8)(a′,e′,ω′,pa′,pe′,pω′)→(a′′,ϕ,ω′′,pa′′,pϕ,pω′′),defined by the following equations: (4.9)a′=a′′,pa′=pa′′,e′=sinϕ,pe′=pϕcosϕ,ω′=ω′′,pω′=pω′′.The Hamiltonian function F′ is invariant with respect to this
transformation. Thus, (4.10)F′′=a′′2μ{4a′′2pa′′2+52pϕ′′2+(52csc2ϕ−2)pω′′2}.The canonical system governed by the new
Hamiltonian functionF′′ has
three first integrals (t and ω are
ignorable variables), (4.11)pω′′=C1,(4.12)pϕ2+pω′′2csc2ϕ=C22,(4.13)a′′2μ{4a′′2pa′′2+52pϕ′′2+(52csc2ϕ−2)pω′′2}=E.(Constant E should
not be confused with the eccentric anomaly E.)
These first integrals play important role in determining the general solution
of the canonical system governed by the new Hamiltonian function F′′ through
Hamilton-Jacobi theory.Let us consider the problem of determining a
canonical transformation defined by a generating functionW such that C1,C2, and E become
the new generalized coordinates, (4.14)(a′′,ϕ,ω′′,pa′′,pϕ,pω′′)→W(C1,C2,E,pC1,pC2,pE).Since the new Hamiltonian function F′′ is a
quadratic form in the conjugate momenta, the separation of variables technique
will be applied for solving the Hamilton-Jacobi equation [3].Consider
the transformation equations(4.15)pa′′=∂W∂a′′,pϕ=∂W∂ϕ,pω′′=∂W∂ω′′,pC1=−∂W∂C1,pC2=−∂W∂C2,pE=−∂W∂E,     where W=W(a′′,ϕ,ω′′,C1,C2,E) is the
generating function. The corresponding Hamilton-Jacobi equation is then given
by (4.16)a′′2μ{4a′′2(∂W∂a′′)2+52(∂W∂ϕ)2+(52csc2ϕ−2)(∂W∂ω′′)2}=E.Following the separation of variables technique,
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, (4.17)W(a′′,ϕ,ω′′,C1,C2,E)=W1(a′′,C1,C2,E)+W2(ϕ,C1,C2,E)+W3(ω′′,C1,C2,E).Therefore,
from (4.11) through (4.17), it follows that (4.18)∂W3∂ω′′=C1,(∂W2∂ϕ)2+(∂W3∂ω′′)2csc2ϕ=C22,a′′2μ{4a′′2(∂W1∂a′′)2+52(∂W2∂ϕ)2+(52csc2ϕ−2)(∂W3∂ω′′)2}=E.A complete solution of these equations is given
by (4.19)W1=−5C22{[4μE5C2a′′−1]1/2−tan−1[4μE5C2a′′−1]1/2},W2=∫C22−C12csc2ϕdϕ,W3=C1ω′′,
with 5C22−4C12=5C2. We note that W2 is given
as indefinite integral, since only its partial derivatives are needed (4.15), as
shown in the next paragraphs.Now,
consider the differential equations for the conjugate momenta of the canonical
system governed by the new HamiltonianF′′′=E, that is, (4.20)dpC1dt=0,dpC2dt=0,dpEdt=−1,  whose solution is very simple: (4.21)pC1=α1,pC2=α2,pE=α3−t,    where αi, i=1,2,3, are arbitrary constants of integration.Introducing
the generating functionW, defined by
(4.17) and (4.19), into the transformation equation (4.15), and taking into
account the general solution defined by (4.21) for the conjugate momenta, one
gets (4.22)pa′′=122a2[4μEa−5C2a2]1/2,pϕ=C22−C12csc2ϕ,pω′′=C1,α1=−ω−225C1Ctan−1((4μE5C2a−1)1/2)+C1∫csc2ϕC22−C12csc2ϕdϕ,α2=−52C2Ctan−1((4μE5C2a−1)1/2)+sin−1(C2C22−C12cosϕ),t−α3=−122Ea[4μEa−5C2a2]1/2.From
equations above, one finds the solution of the canonical system governed by the
HamiltonianF′′ for a given set of initial conditions: (4.23)a′′(t)=a0′′1+(4a0′′/μ)((1/2)Et2−a0′′pa0′′t),(4.24)a′′sin2k0=a0′′sin2(2ψ+k0),(4.25)ψ=15(τ−τ0)1+4cos2k1,(4.26)cosϕ=cosk1cosτ,(4.27)ω′′=k2+tan−1(tanτcsck1)−45τsink1,(4.28)pa′′2=(a0′′a′′)3pa0′′2+18pω0′′2(5csc2k1−4)(a0′′a′′3−1a′′2),(4.29)pϕ2=pω0′′2(csc2k1−csc2ϕ),(4.30)pω′′=pω0′′, with the auxiliary constants k0, k1, and k2 defined
as functions of the initial value of the adjoint variables (conjugate momenta)
by (4.31)csc2k0=8(a0′′pa0′′)2+pω0′′2(5csc2k1−4)pω0′′2(5csc2k1−4),csc2k1=pϕ02+pω0′′2csc2ϕ0pω0′′2,k2=ω0′′+45τ0sink1−tan−1(tanτ0csck1).The constants C, C1, C2, and E in (4.22)  can also be written as
functions of the initial value of the adjoint variables as follows: (4.32)C2=15pω0′′2(5csc2k1−4),C1=pω0′′,C22=pϕ02+pω0′′2csc2ϕ0,4μE=a0′′(8(a0′′pa0′′)2+pω0′′2(5csc2k1−4)). The
initial conditions for the state variables (generalized coordinates) are given
by a′′(0)=a0′′, e′′(0)=sinϕ0, and ω′′(0)=ω0′′,
and τ0 is obtained from cosϕ0=cosk1cosτ0.Equation (4.23) through (4.30) represents
the solution of the canonical system concerning the problem of optimal long
duration low-thrust limited-power transfers between arbitrary elliptic coplanar
orbits in a Newtonian central field.For maneuvers with arbitrary duration, the
periodic terms must be included [9]. Following Hori method and applying the
initial conditions, one gets a first-order analytical solution: (4.33)a(t)=a′(t)+a′5μ3[8e′sinE′a′2p′a+4(1−e′2)sinE′a′pe′−4(1−e′2)1/2e′cosE′a′pω′]E0′E′e(t)=e′(t)+a′5μ3[4(1−e′2)sinE′a′pa′+(1−e′2)[−54e′sinE′+34sin2E′−112e′sin3E′]pe′+(1−e′2)1/2e′[54e′cosE′+14(e′2−3)cos2E′+112e′cos3E′]pω′]E0′E′,ω(t)=ω′(t)+a′5μ3[−4(1−e′2)1/2e′cosE′a′pa′+(1−e′2)1/2e′[54e′cosE′+14(e′2−3)cos2E′+112e′cos3E′]pe′+1e′2[(54−e′2)e′sinE′−12(32−e′2)sin2E′+112e′sin3E′]pω′]E0′E′,
with a′,e′,…,pω′ given by
(4.9). As mentioned before, these equations cannot be applied to transfers
involving very small eccentricities.A
complete first-order analytical solution for the consumption variable is
obtained through a different procedure as described in [9]. This procedure is
based on the linearized solution for transfers between close orbits and it
involves the generating function S.
The expression of J can be put in a
compact form as follows: (4.34)J=E(t−t0)+ΔS,with ΔS=S(E′)−S(E0′).Equations
(4.33) are  in
agreement with the ones obtained in [9] through Bogoliubov-Mitropolsky method, a
noncanonical perturbation technique. But it should be noted that in [9] terms in cos2E′ appear instead of terms in cos2E′;
the other terms are equivalent. The eccentric anomaly E′ is computed from Kepler's equation with the
mean anomaly M′ given by (4.35)M′(t)=M0′+∫t0t[μa′3−(5+2e′22)a′5μ31−e′2e′2pω′]dt,where M0′=M′(t0). The differential equation for the mean
anomaly M′ is
derived from the undisturbed Hamiltonian F0 and an additional term ΔF′ factored by pM′, (4.36)F′=(n′−(5+2e′22)a′5μ31−e′2e′2pω′)pM′+⋯,where dots denote the terms defined by (4.6).For
transfers between coaxial orbits, with no change in the pericenter argument,
the analytical solution described in the preceding paragraphs simplifies. Sincepω0′=0,
some equations must be rewritten: (4.37)a′′sin2k0=a0′′sin2(25(ϕ−ϕ0)+k0),(4.38)pa′′2=(a0′′a′′)3pa0′′2+58pϕ2(a0′′a′′3−1a′′2),(4.39)pϕ=pϕ0=C, with csc2k0=(8(a0′′pa0′′)2+5pϕ2)/5pϕ2 and 4μE=a0′′(8(a0′′pa0′′)2+5pϕ2).Figures3, 4, and 5 show a comparison between the complete first-order analytical
solution, the average solution, and a numerical solution. The numerical
solution has been obtained through integration of the set of canonical
equations which describes the optimal trajectories (system of differential
equations governed by new Hamiltonian function H*, defined by (3.21)–(3.22)). Two sets
of initial conditions and transfer duration defined in Table 1 (in canonical
units) are used in the comparison. In Table 2, final values of state variables
are shown. 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 [10, 11], has been used in the numerical integrations. Note
that there exists an excellent agreement between the complete analytical
solution and the numerical one. Figures 3 and 4 show that the difference
between these solutions is lesser than 10−4 for semimajor axis and eccentricity; but the
analytical solution for pericenter argument loses accuracy for large values of
time. Figure 5 shows the results for the consumption variable. On the other
hand, by comparing the analytical or numerical solution with the average
solution, we see that the amplitudes of the short periodic terms are small, but
they can be significant for transfers with moderate duration. Note that the
difference between the numerical and average solutions decrease as the transfer
duration increases. The good agreement between the numerical and the complete
analytical solution suggest that the latter can be used in the solution  of the
two-point  boundary value problem of going from an initial orbit to a final
orbit, as described in Section 6.Table 1
Set of initial conditions and transfer
duration.
Set 1Set 2a01.01.0e00.200.20ω00.00.0pa02.9326×10−41.4663×10−4pe02.9625×10−51.4812×10−5pω00.00.0tf−t05001000Table 2
Final state variables.
Final  state variablesAnalytical theoryNumerical solutionSecular solutionComplete solutiontf−t0500100050010005001000af2.00002.00001.99671.99731.99681.9973ef0.25000.25000.24750.24800.24760.2480ωf0.00000.0000−0.0094−0.0010−0.0237−0.0081Jf8.5417×10−54.2913×10−58.6327×10−54.3182×10−58.6332×10−54.3182×10−5Comparison between numerical, analytical, and
average solutions for a transfer between coaxial orbits—set 1 of initial
conditions andtf−t0=500.
(a)(b)(c)(d)(e)(f)Comparison between numerical, analytical, and
average solutions for a transfer between coaxial orbits—set 2 of initial
conditions andtf−t0=1000.
(a)(b)(c)(d)(e)(f)Comparison between numerical, analytical, and
average solutions for a transfer between coaxial orbits—consumption variable.
(a)(b)(c)(d)(e)(f)


## 5. Analysis of Long Duration Transfers

In
this section, a complete analysis of long duration transfers described by the
HamiltonianF′ is presented. The investigation of conjugate
points and the solution of the two-point boundary value problem are considered.To
perform this analysis, a new consumption variable is introduced [4], u=2Et,
and (4.23) and (4.24) can be put in the form (5.1)a′′a0′′=[1−2(uv0)cosk0+(uv0)2]−1,(5.2)(uv0)=sin(2ψ)sin(2ψ+k0),where v0=μ/a0′′. By eliminating k0 from
these equations, one finds (5.3)(uv0)=1−2a0′′a′′cos(2ψ)+(a0′′a′′).On the other hand, J=Et, thus (5.4)J=v022t[1−2a0′′a′′cos(2ψ)+(a0′′a′′)].Note that t0=0 in equations above.Now,
consider (4.24)–(4.27) which define
a two-parameter family of extremals in the phase space (a′′,ϕ,ω′′) for a given initial phase point (a0′′,ϕ0,ω0′′) corresponding to an initial orbit. By
eliminating the auxiliary variables τ and ψ, α=a′′/a0′′ and ω′′ can be
written as explicit functions of ϕ, that is, α=α(ϕ,ϕ0,k0,k1) and ω′′=ω′′(ϕ,ϕ0,ω0,k1). The conjugate points to the phase point (a0′′,ϕ0,ω0′′) are determined through the roots of the
equation [12] (5.5)∂α∂k0∂ω′′∂k1−∂α∂k1∂ω′′∂k0=0.Since ω′′ does not
depend on k0, ∂ω′′/∂k0=0. On the other hand, from (4.24) and (5.2), one
finds (5.6)∂α∂k0=−2αsink0(uv0)≠0Thus, the problem of determining the conjugate
points reduces to the analysis of the roots of the following equation: (5.7)∂ω′′∂k1=0.From (4.26) and (4.27), one finds the
explicit form of ω′′=ω′′(ϕ,ϕ0,ω0,k1),
that is, (5.8)ω′′=ω0′′−45sink1{cos−1(cosϕcosk1)−cos−1(cosϕ0cosk1)}+cos−1(cotϕcotk1)−cos−1(cotϕ0cotk1).The partial derivative ∂ω′′/∂k1 is given
by (5.9)∂ω′′∂k1=−45cosk1(τ−τ0)−45tan2k1[cosϕ0sinτ0−cosϕsinτ]+sec2k1[1tan2ϕ0−tan2k1−1tan2ϕ−tan2k1]. The auxiliary variable τ is
reintroduced to simplify the expression.Therefore,
from (5.7) and (5.9), after some simplifications, it follows that the conjugate
points are given by the roots of following equation [4]: (5.10)(1−45sin2k1)sin(τ−τ0)−45cos2k1(τ−τ0)sinτsinτ0=0.As
discussed in [13, 14], the conjugate points determined through (5.10) occur for
transfers during which the direction of motion in the orbit reverses; that is,
for transfers between direct orbits (no direction reversals) the extremals
computed for the canonical system governed by the Hamiltonian F′ are globally minimizing.Some
numerical results are presented through Figures6–9. Figures 6 and 7
show the field of extremals and isocost curves for long duration transfers
between coaxial orbits, considering two cases: in the first one, the space
vehicle departs from a circular orbit and, in the second one, it departs from
an elliptic orbit with eccentricity  e0=0.4. In both cases, the semimajor axis is a0=1,
such that the results are presented in canonical units. The extremals are
plotted for different values of the constant k0 and the isocost curves for different values of
the new consumption variable u/v0.
Only transfers between direct orbits are considered.Figure 6
Field of extremals and isocost curves for
transfers  between coaxial orbits fore0=0.0.Figure 7
Field of extremals and isocost curves for
transfers   between coaxial orbits fore0=0.4.Field of extremals for transfers  between noncoaxial orbits fore0=3/2
with k0=18°.
(a)(b)Figure 9
Isocost curves fora0=1.0, e0=0.5, and ω0=0°,
with k1=18°.Figures8 and 9 show the field of extremals and isocost curves for long duration transfers
between noncoaxial orbits. Note that in the general case the family of
extremals is described by two constants k0 and k1.
In the plots of Figure 8, k0=18° has been chosen in order to illustrate the
general aspect of the field of extremals and compare the results with the ones
obtained in [4, 5, 13, 14]. The projection of extremals on the (e,ω) -plane shows clearly the nonexistence of
conjugate points for the extremals computed from (4.6)–(4.11) in the case
of transfers between direct orbits (no direction reversals).Figure9 shows a sample of isocost curves for different values of the new consumption
variable u/v0 and a fixed value of k1=18°.
Note that these curves lie on an isocost surface with complex geometry. For
simplicity, the symmetric interval [−180°,180°] for ω is considered in this representation.To complete the analysis, a brief
description of an iterative algorithm for solving the two-point boundary value
problem of going from an initial orbitO0:(a0′′,e0′′,ω0′′) to a final orbit Of:(af′′,ef′′,ωf′′) is presented. The steps of this algorithm can
be described as follows.(1)
Guess a starting value ofk1,(2)
Determineϕ0,ϕf,τ0, and τf through (4.9) and (4.26).(3)
Computeω′′ through (4.27).(4)
Ifω′′≠ωf′′,
adjust the value of k1 and repeat steps (2) and (3) until |ω′′−ωf′′|<δ,
where δ is a prescribed small quantity.(5)
Computeψf and (uf/v0) through (4.25) and (5.3).(6)
Computek0 using equation tank0=sin(2ψf)/(αf−cos(2ψf)),
where αf=af′′/a0′′.(7)
Compute  successivelypa0′′, pω0′′, and pϕ0′′ through the  following  equations: (5.11)pa0′′=v022a0′′tf(ufv0)cosk0,pω0′′=22a0′′pa0′′cotk05csc2k1−4,pϕ02=pω0′′2(csc2k1−csc2ϕ0). The
adjust of the value of k1 can be obtained through the well-known
Newton-Raphson method, that is, (5.12)k1n+1=k1n−ω′′(k1n)−ωf′′(∂ω′′/∂k1)k1=k1n, with
the partial derivative (∂ω′′/∂k1)k1=k1n computed from (5.9).In the case of coaxial orbits,k1=0 and the solution of the two-point boundary
value problem involves only the solution of (4.37), which is given through step (6)
with ψ replaced by (1/5)(ϕ−ϕ0).
The initial value of the adjoint variables is given as follows: pa0′′ is computed through step (7) and pϕ0 is given by (5.13)pϕ02=8tan2k0(a0′′pa0′′)25.With the solution of the two-point
boundary value problem, the optimal thrust acceleration required to perform the
maneuver can be determined through (2.5) and (3.22), that is, (5.14)γ*=1na1−e2{{2aesinfpa+((1−e2)sinf)pe−(1−e2)cosfepω}er+{2a(1−e2)(ar)pa+(1−e2)(cosf+cosE)pe+(1−e2)sinfe(1+11+ecosf)pω}es}. The
state and adjoint variables are computed from (4.23) through (4.30). For simplicity,
prime denoting the average variables is omitted.


## 6. Solution of the Two-point Boundary Value Problem

In
this section, an iterative algorithm based on the approximate analytical
solution defined by (4.33) is briefly described for solving the
two-point boundary value problem of going from an initial orbit O0:(a0,e0,ω0) to a
final orbit Of:(af,ef,ωf) at the
prescribed final time tf.For
a given time, (4.33) can be represented as follows: (6.1)yi(t)=gi(t,pa′,pe′,pω′),i=1,2,3,where y1(t)=a(t), y2(t)=e(t), and y3(t)=ω(t).
Note that pa′, pe′, and pω′ appear explicitly
in the short periodic terms, but they also appear implicitly through a′, e′, ω′, and E′.
Thus, functions gi, i=1,2,3,
are nonlinear in these variables.So,
the two-point boundary value problem can be stated as follows: findpa′, pe′, and pω′ such
that y1(tf)=af, y2(tf)=ef, and y3(tf)=ωf.
This problem can be solved through a Newton-Raphson algorithm [10, 11] which
uses the algorithm described in the preceding section (for long duration
transfers) to generate a starting approximation of the adjoint variables.Figure10 shows a comparison between the solution of the two-point boundary value
problem based on the average and complete solutions for a transfer between
coaxial orbits with the departure orbit O0 defined by the set of initial conditions a0=1.0, e0=0.20, ω0=0° and the
arrival orbit defined by the set of terminal conditions af=2.0, ef=0.25, ωf=0°, and two transfer durations tf−t0=500 and tf−t0=1000.
All variables are represented in canonical units.Complete and average solutions of the
two-point boundary value problem  fortf−t0=500 and tf−t0=1000.
(a)(b)From
results above, we note that average solution of the TPBVP does not represents a
mean solution of the complete solution of the same TPBVP; that is, these
solutions have different values for the initial adjoint variables, and, the
solution of the TPBVP based on the average solution becomes closer to the
solution based on the complete analytical one, as the transfer duration
increases.


## 7. Transfers between Close Elliptical Coplanar Orbits

For
transfers between close elliptical coplanar orbits, an approximate solution of
the system of differential equations governed by the HamiltonianH∗,
defined through (3.21) and (3.22), can be straightforwardly obtained by linearizing
the functions F1 and S1 about a reference elliptic orbit O¯ defined by a nominal set of orbital elements a¯, e¯, and ω¯.
This solution is very similar to the complete solution given by (4.33), with a¯, e¯, and ω¯ replacing a′, e′, ω′,
and can be put in a suitable matrix form (7.1)Δx=Ap0,where Δx=[ΔαΔeΔω]T denotes the imposed changes on orbital
elements (state variable), α=a/a¯, pα=a¯pa, p0 is the 3×1 matrix of initial values of the adjoint
variable and A is 3×3 symmetric matrix. The adjoint variables are
constant, and the matrix A is given by (7.2)A=[aααaαeaαωaeαaeeaeωaωαaωeaωω], where (7.3)aαα=4a¯5μ3(E¯+e¯sinE¯)|E¯0E¯f,aαe=aeα=4a¯5μ3(1−e¯2)sinE¯|E¯0E¯f,aαω=aωα=−4a¯5μ3(1−e¯2)1/2e¯cosE¯|E¯0E¯f,aee=a¯5μ3(1−e¯2)[52M¯−54e¯sinE¯+34sin2E¯−112e¯sin3E¯]|E¯0E¯f,aeω=aωe=a¯5μ3(1−e¯2)1/2e¯[54e¯cosE¯+14(e¯2−3)cos2E¯+112e¯cos3E¯]|E¯0E¯f,aωω=a¯5μ31e¯2[(52−2e¯2)M¯+e¯(54−e¯2)sinE¯−12(32−e¯2)sin2E¯+112e¯sin3E¯]|E¯0E¯f,
where E¯ is calculated by solving Kepler’s equation for M¯=n¯(t−τ),
with n¯=μ/a¯3.
Here, τ is the time of pericenter passage. The
overbar denotes the reference orbit.Equations
(7.1)–(7.3) represent an
approximate first-order solution for optimal low-thrust limited-power transfers
between close elliptic coplanar orbits. The initial value of the adjoint
variables must be determined to satisfy the two-point boundary value problem of
going from an initial orbit O0:(a0,e0,ω0) to a
final orbit Of:(af,ef,ωf) at the
prescribed final time tf.
In this case, the vector p0 can be determined by simple techniques since
(7.1) defines a linear system of algebraic equations [15–17]; no iterative
method such as Newton-Raphson is needed.Figure11 shows a comparison between the complete (nonlinear) first-order solution
defined by (4.33) and the
approximate (linearized) first-order solution defined by (7.1)–(7.3) for three
set of terminal orbits described in Table 3, with tf−t0=5 and
100 canonical units. Note that linearized solution provides a good
approximation for the nonlinear first-order solution. This approximation is
better for transfers with larger duration. For transfers with moderate duration
corresponding to some revolutions, the amplitude of periodic terms decreases as
the eccentricities of terminal orbits increase.Table 3
Set of terminal orbits.
Initial orbit 1Final orbit 1Initial orbit 2Final orbit 2Initial orbit 3Final orbit 3a1.0001.0501.0001.0501.0001.050e0.1000.1050.2000.2100.5000.525ω0.00.00.00.00.00.0Linearized and nonlinear solutions for
transfers between close orbits.
(a)(b)(c)(d)(e)(f)


## 8. Conclusion

A
first-order analytical solution for the problem of optimal low-thrust
limited-power transfers between arbitrary elliptic coplanar orbits in a Newtonian
central gravity field has been obtained through Hamilton-Jacobi theory and a
perturbation method based on Lie series. A complete analysis of long duration
transfers, which includes the investigation of conjugate points and the
solution of the two-point boundary value problem, is presented. Results show
the good agreement between the complete analytical solution and the one
obtained through numerical integration of the set of canonical equations which
describes the optimal trajectories for a transfer between coaxial orbits. On
the other hand, the average solution provides a good approximation for large
values of the time when compared to the complete analytical solution. For
transfers between close coplanar elliptic orbits, a linearized solution is
straightforwardly obtained by linearizing the new Hamiltonian functionF1 and the generating function S1 about a reference elliptic orbit O¯.
Numerical results show the good agreement between the linearized solution and
the complete (nonlinear) analytical solution for such transfers. Finally, we
note that a first-order solution for the Cartesian elements can also be easily
obtained from Lie theorem [2], since the theory presented in the paper is based
on canonical transformations.



---

*Source: 525930-2008-10-15.xml*