Optimization of Low-Thrust Limited-Power Trajectories in a Noncentral Gravity Field—Transfers between Orbits with Small Eccentricities

Abstract

Numerical and first-order analytical results are presented for optimal low-thrust limited-power trajectories in a gravity field that includes the second zonal harmonic 𝐽₂ in the gravitational potential. Only transfers between orbits with small eccentricities are considered. The optimization problem is formulated as a Mayer problem of optimal control with Cartesian elements—position and velocity vectors—as state variables. After applying the Pontryagin Maximum Principle, successive canonical transformations are performed and a suitable set of orbital elements is introduced. Hori method—a perturbation technique based on Lie series—is applied in solving the canonical system of differential equations that governs the optimal trajectories. First-order analytical solutions are presented for transfers between close orbits, and a numerical solution is obtained for transfers between arbitrary orbits by solving the two-point boundary value problem described by averaged maximum Hamiltonian, expressed in nonsingular elements, through a shooting method. A comparison between analytical and numerical results is presented for some maneuvers.

1. Introduction

This paper presents numerical and analytical results for optimal low-thrust limited-power transfers between orbits with small eccentricities in a noncentral gravity field that includes the second zonal harmonic 𝐽₂ in the gravitational potential. 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 two space mission powered by ionic propulsion: Deep Space One and SMART 1. Several researchers have obtained numerical and sometimes analytical solutions for a number of specific initial orbits and specific thrust profiles in central or noncentral gravity field [1–26].

It is assumed that the thrust direction is free and the thrust magnitude is unbounded [5, 27], that is, there exist no constraints on control variables. It is supposed that the orbital changes caused by the thrust and the perturbations due to the oblateness of the central body have the same order of magnitude. The optimization problem is formulated as a Mayer problem of optimal control with position and velocity vectors—Cartesian elements—as state variables. After applying the Pontryagin Maximum Principle [28], successive Mathieu transformations are performed and a suitable set of orbital elements is introduced. Hori method—a perturbation technique based on Lie series [29]—is applied in solving the canonical system of differential equations that governs the optimal trajectories. For transfers between close quasicircular orbits, new set of nonsingular orbital elements is introduced and the new Hamiltonian and the generating function, determined through the algorithm of Hori method, are transformed to these new elements through Mathieu transformations. First-order analytical solutions are obtained for transfers between close quasicircular orbits, considering maneuvers between near-equatorial orbits or nonequatorial orbits. These solutions are given in terms of systems of algebraic equations involving the imposed variations of the orbital elements, the duration of the maneuvers, the effects of the oblateness of the central body, and the initial values of the adjoint variables. 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 [5]. For transfers between arbitrary orbits, a slightly different set of nonsingular elements is introduced, and the two-point boundary value problem described by averaged maximum Hamiltonian, expressed in nonsingular elements, is solved through a shooting method. A comparison between analytical and numerical results is presented for some maneuvers.

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 [5]. The ratio between the maximum thrust acceleration and the gravity acceleration on the ground, 𝛾_max/𝑔₀, is between 10⁻⁴ and 10⁻². For such system, the fuel consumption is described by the variable 𝐽 defined as 𝐽=12^𝑡𝑡₀𝛾²𝑑𝑡,(2.1) where γ is the magnitude of the thrust acceleration vector 𝜸, used as control variable. The consumption variable 𝐽 is a monotonic decreasing function of the mass 𝑚 of space vehicle: 𝐽=𝑃_max𝑚1−𝑚1₀,(2.2) where 𝑃_max is the maximum power, and 𝑚₀ is the initial mass. The minimization of the final value of the fuel consumption 𝐽_𝑓 is equivalent to the maximization of 𝑚_𝑓.

The general optimization problem concerned with low-thrust limited-power transfers (no rendezvous) will be formulated as a Mayer problem of optimal control by using Cartesian elements as state variables. Consider the motion of a space vehicle 𝑀 powered by a limited-power engine in a general gravity field. At time 𝑡, the state of the vehicle is defined by the position vector 𝐫(𝑡), the velocity vector 𝐯(𝑡), and the consumption variable 𝐽. The geometry of the transfer problem is illustrated in Figure 1. The control 𝜸 is unconstrained, that is, the thrust direction is free, and the thrust magnitude is unbounded.

Figure 1 

Geometry of the transfer problem.

The optimization problem is formulated as follows: it is proposed to transfer the space vehicle 𝑀 from the initial state (𝐫₀,𝐯₀,0) at the initial time 𝑡₀=0 to the final state (𝐫_𝑓,𝐯_𝑓,𝐽_𝑓) at the specified final time 𝑡_𝑓, such that the final consumption variable 𝐽_𝑓 is a minimum. The state equations are 𝑑𝐫𝑑𝑡=𝐯,𝑑𝑑𝐯𝑡=𝐠(𝐫)+𝜸,𝑑𝑑𝐽𝑡=12𝛾²,(2.3) where 𝐠(𝐫) is the gravity acceleration. For transfers between orbits with small eccentricities, the initial and final conditions will be specified in terms of singular orbital elements introduced in next sections. It is also assumed that 𝑡_𝑓−𝑡₀ and the position of the vehicle in the initial orbit are specified.

According to the Pontryagin Maximum Principle [28], the optimal thrust acceleration 𝜸^∗ must be selected from the admissible controls such that the Hamiltonian function 𝐻 reaches its maximum. The Hamiltonian function is formed using (2.3), 𝐻=𝐩_𝑟•𝐯+𝐩_𝑣•(𝐠(𝐫)+𝜸)+12𝑝^𝐽𝛾²,(2.4) where 𝐩_𝑟,𝐩_𝑣, and 𝑝_𝐽 are the adjoint variables, and dot denotes the dot product. Since the optimization problem is unconstrained, 𝜸^∗ is given by 𝜸^∗=−𝐩^𝑣𝑝_𝐽.(2.5) The optimal thrust acceleration 𝜸^∗ is modulated [5], and the optimal trajectories are governed by the maximum Hamiltonian function 𝐻^∗, obtained from (2.4) and (2.5): 𝐻^∗=𝐩_𝑟•𝐯+𝐩_𝑣•𝐠(𝐫)−𝑝^𝑣22𝑝_𝐽.(2.6)

The consumption variable 𝐽 is ignorable, and 𝑝_𝐽 is a first integral. From the transversality conditions, 𝑝_𝐽(𝑡_𝑓)=−1; thus, 𝑝_𝐽(𝑡)=−1.(2.7) Equation (2.6) reduces to 𝐻=𝐩_𝑟•𝐯+𝐩_𝑣•𝐠(𝐫)+𝑝^𝑣22.(2.8) The consumption variable 𝐽 is determined by simple quadrature.

For noncentral field, the gravity acceleration is given by [30] 𝐠(𝐫)=−𝜇𝑟³𝐫+𝜕𝑉²𝜕𝐫^𝑇,(2.9) with the disturbing function 𝑉₂ defined by 𝑉₂=−𝜇𝐽₂𝑎𝑟^𝑒23𝑃²(sin𝜑),(2.10) where 𝜇 is the gravitational parameter, 𝑎_𝑒 is the mean equatorial radius of the central body, 𝐽₂ is the coefficient for the second zonal harmonic of the potential, 𝜑is the latitude, and 𝑃₂ is the Legendre polynomial of order 2. (𝜕𝑉₂/𝜕𝐫) is a row vector of partial derivatives, and superscript 𝑇 denotes its transpose.

Using (2.8) and (2.9), the maximum Hamiltonian function can be put in the following form: 𝐻^∗=𝐻₀+𝐻_𝐽2+𝐻_𝛾∗,(2.11) where 𝐻₀=𝐩_𝑟•𝐯−𝐩_𝑣•𝜇𝑟³𝐻𝐫,_𝐽2=𝐩_𝑣•𝜕𝑉²𝜕𝐫^𝑇𝐻_𝛾∗,=𝑝^𝑣22.(2.12)𝐻₀ is the undisturbed Hamiltonian function, 𝐻_𝐽2 and 𝐻_𝛾∗ are disturbing functions concerning the oblateness of the central body and the optimal thrust acceleration, respectively. 𝐻_𝐽2 and 𝐻_𝛾∗ are assumed to be of the same order of magnitude.

3. Transformation from Cartesian Elements to a Set of Orbital Elements

Consider the canonical system of differential equations governed by the undisturbed Hamiltonian 𝐻₀: 𝑑𝐫𝑑𝑡=𝐯,𝑑𝑑𝐯𝑡=−𝑟𝜇³𝑑𝐩_𝐫𝐫,𝑑𝑡=𝑟𝜇³𝐩_𝑣−3𝐩_𝑣•𝐞_𝐫𝐞_𝐫,𝑑𝐩^𝑣𝑑𝑡=−𝐩^𝑟,(3.1) where 𝐞_𝐫 is the unit vector pointing radially outward of the moving frame of reference (Figure 2). The general solution of the state equations is well known in Astrodynamics [31], and the general solution of the adjoint equations is obtained through properties of generalized canonical systems [32], as described in the appendix. Thus, 𝐫=𝑎1−𝑒²1+𝑒cos𝑓𝐞^𝐫𝐯=,(3.2)𝜇𝑎1−𝑒²𝑒sin𝑓𝐞^𝐫𝐩+(1+𝑒cos𝑓)𝐞_𝐬,(3.3)_𝐫=𝑎𝑟²2𝑎𝑝^𝑎+1−𝑒²cos𝐸𝑝_𝑒+𝑟𝑎sin𝑒𝑓𝑝^𝜔−1−𝑒³cos𝐸√1−𝑒²𝑝_𝑀𝐞^𝐫+sin𝑓𝑎𝑝^𝑒−(𝑒+cos𝑓)𝑎𝑒1−𝑒²𝑝^𝜔+√1−𝑒²cos𝑓𝑎𝑒𝑝^𝑀𝐞^𝐬+1𝑎√1−𝑒²+√𝑎𝑟sin𝐸𝑝^𝐼cos𝜔+si𝑝n^Ω𝐼−𝑝^𝜔cot𝐼sin𝜔1−𝑒²𝑎𝑟cos𝐸𝑝^𝐼sin𝜔−𝑝^Ωsin𝐼−𝑝^𝜔𝐩cot𝐼cos𝜔𝐞_𝐰,(3.4)_𝐯=1𝑛𝑎√1−𝑒²2𝑎𝑒sin𝑓𝑝^𝑎+1−𝑒²sin𝑓𝑝_𝑒−1−𝑒²cos𝑓𝑒𝑝^𝜔+1−𝑒²^3/2𝑒cos𝑓−1+𝑒2𝑒cos𝑓𝑝^𝑀𝐞^𝐫+2𝑎1−𝑒²𝑎𝑟𝑝^𝑎+1−𝑒²(cos𝑓+cos𝐸)𝑝_𝑒+1−𝑒²sin𝑓𝑒1+1+𝑒1cos𝑓𝑝^𝜔−√1−𝑒²𝑝_𝑀𝐞^𝐬+𝑟𝑎cos(𝜔+𝑓)𝑝^𝐼+𝑟𝑎sin(𝜔+𝑓)si𝑝n^Ω𝐼−𝑝^𝜔cot𝐼𝐞_𝐰,(3.5) where 𝐞_𝐬 and 𝐞_𝐰 are unit vectors along circumferential and normal directions of the moving frame of reference, respectively; 𝑎 is the semimajor axis, 𝑒 is the eccentricity, 𝐼 is the inclination of orbital plane, Ω is the longitude of the ascending node, 𝜔is the argument of pericenter, 𝑓 is the true anomaly, 𝐸 is the eccentric anomaly, 𝑀 is the mean anomaly, 𝑛=√𝜇/𝑎³ is the mean motion, and (𝑟/𝑎),(𝑟/𝑎)sin𝑓, and so forth are functions of the elliptic motion which can be expressed explicitly in terms of the eccentricity and the mean anomaly through Lagrange series [31]. The true, eccentric and mean anomalies are related through the equations: tan𝑓2=1+𝑒𝑀=𝐸−1𝑒−s𝑒inta𝐸n.𝐸2,(3.6) The unit vectors 𝐞_𝐫,𝐞_𝐬, and 𝐞_𝐰 of the moving frame of reference are written in the fixed frame of reference as 𝐞𝐞_𝐫=+(c(ossinΩΩcocsos(𝜔(𝜔++𝑓𝑓))−+sicnoΩsΩsinsin(𝜔(𝜔++𝑓)𝑓c)ocso𝐼s)𝐼𝐢)𝐣+sin(𝜔+𝑓)sin𝐼𝐤,_𝐬𝐞=+−((−cossinΩΩsisnin(𝜔(𝜔++𝑓𝑓))++sicnoΩsΩcocsos(𝜔(𝜔++𝑓)𝑓c)ocso𝐼s)𝐼𝐢)𝐣+cos(𝜔+𝑓)sin𝐼𝐤,_𝐰=sinΩsin𝐼𝐢−cosΩsin𝐼𝐣+cos𝐼𝐤.(3.7)

Figure 2 

Frames of reference.

Equations (3.2)–(3.5) define a Mathieu transformation between the Cartesian elements (𝐫,𝐯,𝐩_𝐫,𝐩_𝐯) and the orbital ones (𝑎,𝑒,𝐼,Ω,𝜔,𝑀,𝑝_𝑎,𝑝_𝑒,𝑝_𝐼,𝑝_Ω,𝑝_𝜔,𝑝_𝑀). Since the Hamiltonian function is invariant with respect to this canonical transformation, it follows, from (2.11) through (3.5), that 𝐻^∗=𝐻₀+𝐻_𝐽2+𝐻_𝛾∗,(3.8) where 𝐻𝐻_0𝐽2=𝑛𝑝_𝑀,=2𝑛𝑎𝜕𝜕𝑀𝑉^2𝑝^𝑎+√1−𝑒²𝑛𝑎²𝑒−𝜕𝜕𝑉𝜔^2+√1−𝑒²𝜕𝑉^2𝜕𝑀𝑝^𝑒+1𝑛𝑎²√1−𝑒²sin𝐼−𝜕𝜕𝑉Ω^2+cos𝐼𝜕𝜕𝑉𝜔^2𝑝^𝐼+1𝑛𝑎²√1−𝑒²sin𝐼𝜕𝜕𝑉𝐼^2𝑝^Ω+√1−𝑒²𝑛𝑎²𝑒𝜕𝜕𝑉𝑒^2−𝑒1c−ot𝑒𝐼²𝜕𝜕𝑉𝐼^2𝑝^𝜔+1𝑛𝑎−2𝜕𝜕𝑉𝑎^2−1−𝑒²𝐻_∗𝛾𝑎𝑒𝜕𝜕𝑉𝑒^2𝑝^M,=12𝑛²𝑎²1−𝑒²21(1−cos2𝑓)2𝑎𝑒𝑝^𝑎+1−𝑒²𝑝_𝑒²+21−𝑒²sin2𝑓−𝑎𝑝_𝑎𝑝_𝜔−1−𝑒²2𝑒𝑝^𝑒𝑝_𝜔+41−𝑒₂^3/2sin𝑓1+−𝑒2c𝑒os𝑓+cos𝑓𝑎𝑝^𝑎𝑝_𝑀+1−𝑒²2𝑒𝑝^𝑒𝑝_𝑀+1−𝑒₂²2𝑒²(1+cos2𝑓)𝑝_𝜔2−21−𝑒²^5/2𝑒²1+−𝑒2c𝑒os𝑓+cos𝑓cos𝑓𝑝^𝜔𝑝_𝑀+1−𝑒²³𝑒²1+−𝑒2c𝑒os𝑓+cos𝑓²𝑝_𝑀2+4𝑎²1−𝑒²²𝑎𝑟²𝑝_𝑎2+4𝑎1−𝑒²²𝑎𝑟(cos𝐸+cos𝑓)𝑝^𝑎𝑝_𝑒+1−𝑒²²(cos𝐸+cos𝑓)²𝑝_𝑒2+4𝑎1−𝑒²²×𝑝𝑒𝑎𝑟sin𝑓1+1+𝑒1cos𝑓_𝑎+21𝑝−^𝜔−𝑒₂1−𝑒²^1/2𝑝_𝑎𝑝_𝑀²×sin𝑓𝑒𝑝_𝑒(cos𝐸+cos𝑓)1+1+𝑒1cos𝑓+1−𝑒₂𝑝_𝜔−√1−𝑒²𝑝_𝑒𝑝_𝑀𝑒1+1+𝑒1cos𝑓sin𝑓𝑝^𝜔−√1−𝑒²𝑝_𝑀²+12𝑎𝑟²𝑝^𝐼2+𝑝^Ωsin𝐼−𝑝^𝜔cot𝐼²+12𝑎𝑟²cos2(𝜔+𝑓)𝑝_𝐼2−𝑝^Ωsin𝐼−𝑝^𝜔cot𝐼²+𝑟𝑎²sin2(𝜔+𝑓)𝑝_𝐼𝑝^Ωsin𝐼−𝑝^𝜔cot𝐼,(3.9) with the disturbing function 𝑉^2 given by 𝑉^2=𝜇𝑎𝑎𝑎^𝑒²𝐽₂𝑎𝑟³12−34sin²𝐼+34𝑎𝑟³sin²𝐼cos2(𝜔+𝑓).(3.10)

The new Hamiltonian function 𝐻^∗, defined by (3.8)–(3.10), describes the optimal low-thrust limited-power trajectories in a noncentral gravity field which includes the perturbative effects of the second zonal harmonic of the gravitational potential.

4. Averaged Maximum Hamiltonian for Optimal Transfers

In order to eliminate the short periodic terms from the maximum Hamiltonian function 𝐻^∗, Hori method [29] is applied. It is assumed that 𝐻₀ is of zero-order and 𝐻^∗𝛾 is of the first-order in a small parameter associated to the magnitude of the thrust acceleration, and 𝐻_𝐽2 has the same order of 𝐻^∗𝛾.

Consider an infinitesimal canonical transformation: 𝑎,𝑒,𝐼,Ω,𝜔,𝑀,𝑝_𝑎,𝑝_𝑒,𝑝_𝐼,𝑝_Ω,𝑝_𝜔,𝑝_𝑀⟶𝑎^,𝑒^,𝐼^,Ω^,𝜔^,𝑀^,𝑝^𝑎,𝑝^𝑒,𝑝^𝐼,𝑝^Ω,𝑝^𝜔,𝑝^𝑀.(4.1) The new variables are designated by the prime. According to the algorithm of Hori method [29], at order 0, 𝐹₀=𝑛^𝑝^𝑀.(4.2)𝐹₀ denotes the new undisturbed Hamiltonian. Now, consider the canonical system described by 𝐹₀: 𝑑𝑎^𝑑Ω_𝑑𝑡=0,𝑑𝑑𝑒𝑡^=0,𝑑𝑑𝐼𝑡^=0,𝑑𝑡=0,𝑑𝑑𝜔𝑡^=0,𝑑𝑑𝑀𝑡^=𝑛^𝑑𝑝_𝑎,𝑑𝑡=32𝑎𝑛^𝑝^𝑀,𝑑𝑝^𝑒𝑑𝑝_Ω𝑑𝑡=0,𝑑𝑑𝑝𝑡^𝐼=0,𝑑𝑡=0,𝑑𝑑𝑝𝑡^𝜔=0,𝑑𝑑𝑝^𝑀𝑡=0,(4.3) general solution of which is given by Ω_𝑎^=𝑎^0,𝑒^=𝑒^0,𝐼^=𝐼^0,𝑝_𝑎=Ω^0,𝜔^=𝜔^0,𝑀^=𝑀^0+𝑛^𝑡−𝑡₀,=𝑝^𝑎₀+32𝑛^𝑡−𝑡₀𝑎^𝑝_Ω𝑝^𝑀,𝑝^𝑒=𝑝^𝑒₀,𝑝^𝐼=𝑝^𝐼₀,=𝑝^Ω₀,𝑝^𝜔=𝑝^𝜔₀,𝑝^𝑀=𝑝^𝑀₀.(4.4) The subscript 0 denotes the constants of integration.

Introducing the general solution defined above into the equation of order 1 of the algorithm of Hori method, it reduces to 𝜕𝑆₁𝜕𝑡=𝐻^𝐽2+𝐻_𝛾∗−𝐹₁.(4.5) According to [29], the mean value of 𝐻_𝐽2+𝐻_𝛾∗ must be calculated, and, 𝑆₁ is obtained through integration of the remaining part. 𝐹₁ and 𝑆₁ are given by the following equations: 𝐹₁=−32𝑛^𝐽₂𝑎^𝑒𝑎^²1−𝑒^2⁻²cos𝐼^𝑝^Ω+121−5cos²+𝑎_𝐼^𝑝^𝜔2𝜇4𝑎^2𝑝^𝑎2+521−𝑒^2𝑝^𝑒2+5−4𝑒^22𝑒^2𝑝^𝜔2+𝑝^𝐼221−𝑒^21+32𝑒^2+25𝑒^2cos2𝜔^+5𝑒^2sin2𝜔^21−𝑒^2𝑝^𝐼si𝑝n^Ω𝐼^−cot𝐼^𝑝^𝜔+121−𝑒^2si𝑝n^Ω𝐼^−cot𝐼^𝑝^𝜔²1+32𝑒^2−52𝑒^2𝑆₁cos2𝜔^,=𝐽₂𝑎^𝑒𝑎^²1−32sin²𝐼^𝑎^𝑟^³−1−𝑒^2^−3/2+32sin²𝐼^𝑎𝑟^³cos2𝜔^+𝑓^𝑎^𝑝^𝑎+32sin²𝐼^𝑒^1−𝑒^2−12cos2𝜔^+𝑓^−𝑒^2cos2𝜔^+𝑓^+13cos2𝜔^+1−𝑒_2+3𝑓^𝑒^12−34sin²𝐼^𝑎^𝑟^³−1−𝑒^2^−3/2+34sin²𝐼^𝑎𝑟^³cos2𝜔^+𝑓^𝑝^𝑒+34sin2𝐼^1−𝑒^2²21cos2𝜔^+𝑓^+𝑒^2cos2𝜔^+𝑓^+13cos2𝜔^+3𝑓^𝑝^𝐼+32cos𝐼^1−𝑒^2²+𝑒_−𝑓^−𝑀^+𝑒^sin𝑓^+12sin2𝜔^+𝑓^2sin2𝜔^+𝑓^+13sin2𝜔^+3𝑓^𝑝^Ω+345cos²𝐼^−11−𝑒^2²+1𝑓^−𝑀^+𝑒^sin𝑓^43cos²𝐼^−1𝑒^1−𝑒^2𝑎𝑟^²1−𝑒^2+𝑎𝑟^+1sin𝑓^+38sin²𝐼^𝑒^1−𝑒^2−𝑎𝑟^²+𝑎_1−𝑒^2−𝑎𝑟^+1sin2𝜔^+𝑓^𝑟^²1−𝑒^2+𝑎𝑟^+31sin2𝜔^+3+3𝑓^83−5cos²𝐼^1−𝑒^2²+𝑒_sin2𝜔^+𝑓^+1sin2𝜔^+𝑓^+31sin2𝜔^+3𝑓^𝑝^𝜔2𝑎^5𝜇³8𝑒^sin𝐸^𝑎^2𝑝^𝑎2+81−𝑒^2sin𝐸^𝑎^𝑝^𝑎𝑝^𝑒−8√1−𝑒^2𝑒^cos𝐸^𝑝^𝑎𝑝^𝜔+1−𝑒^2−54𝑒^sin𝐸^+34sin2𝐸^−112𝑒^sin3𝐸^𝑝^𝑒2+√1−𝑒^2𝑒^52𝑒^cos𝐸^−123−𝑒^2cos2𝐸^+16𝑒^cos3𝐸^𝑝^𝑒𝑝^𝜔+1−𝑒^2⁻¹𝑝^𝐼2+𝑝^Ωsin𝐼^−𝑝^𝜔cot𝐼^×−𝑒_²+38𝑒^3sin𝐸^+38𝑒^2sin2𝐸^−124𝑒^3sin3𝐸^+1−𝑒^2⁻¹×𝑝^𝐼2cos2𝜔^+2𝑝^𝐼si𝑝n^Ω𝐼^−𝑝^𝜔cot𝐼^sin2𝜔^−𝑝^Ωsin𝐼^−𝑝^𝜔cot𝐼^²cos2𝜔^×−54𝑒^+58𝑒^3sin𝐸^+14+18𝑒^2sin2𝐸^+−112𝑒^+1+1−𝑒_224𝑒^3sin3𝐸^^−1/2−𝑝^𝐼2sin2𝜔^+2𝑝^𝐼𝑝^Ωsin𝐼^−𝑝^𝜔cot𝐼^cos2𝜔^+𝑝^Ωsin𝐼^−𝑝^𝜔cot𝐼^²sin2𝜔^×54𝑒^cos𝐸^−14+14𝑒^2cos2𝐸^+112𝑒^+𝑝_𝜔2cos3𝐸^𝑒^254𝑒^−𝑒^3sin𝐸^+−34+21𝑒^2sin2𝐸^+112𝑒^sin3𝐸^.(4.6) Terms factored by 𝑝^𝑀 have been omitted in equations above, since only transfers (no rendez-vous) are considered [2].

It should be noted that the averaged maximum Hamiltonian and the generating function become singular for circular and/or equatorial orbits. In order to avoid these singularities, suitable sets of nonsingular elements will be introduced in the next sections.

5. Optimal Transfers between Close Quasicircular Orbits

In this section, approximate first-order analytical solutions will be obtained for transfers between close quasicircular orbits.

5.1. Transfers between Close Quasicircular Nonequatorial Orbits

Consider the Mathieu transformation [33] defined by 𝑎^𝐼=_𝑎^,ℎ^=𝑒^cos𝜔^,𝑘^=𝑒^sin𝜔^,𝑝=𝐼^,Ω^=Ω^,ℓ^=𝑀^+𝜔^,(5.1)_𝑎=𝑝^𝑎,𝑝^ℎ=𝑝^𝑒cos𝜔^−𝑝^𝜔−𝑝^𝑀𝑒^sin𝜔^,𝑝^𝑘=𝑝^𝑒sin𝜔^+𝑝^𝜔−𝑝^𝑀𝑒^cos𝜔^𝑝_𝐼,=𝑝^𝐼,𝑝^Ω=𝑝^Ω,𝑝^ℓ=𝑝^𝑀.(5.2) The new set of canonical variables, designated by the double prime, is nonsingular for circular nonequatorial orbits.

Introducing (5.2) into (4.6) and using the expansions of the elliptic motion in terms of eccentricity and mean anomaly [31], one finds that the new averaged maximum Hamiltonian and the generating function are written up to the zeroth order in eccentricity as 𝐹₁=−32𝑛𝐽²𝑎𝑎^𝑒²cos𝐼𝑝_Ω+𝑎2𝜇4𝑎²𝑝^2𝑎+52𝑝^2ℎ+𝑝^2𝑘+12𝑝^2𝐼+𝑝^2Ω𝑆sin²𝐼₁,=32𝐽²𝑎𝑎^𝑒²𝑎𝑝^𝑎sin²𝐼cos2ℓ+𝑝_ℎcosℓ+14sin²𝐼−5cosℓ+7+𝑝_𝑘3cos3ℓ+1sinℓ+41sin²𝐼−7sinℓ+73sin3ℓ4𝑝^𝐼sin2𝐼cos2ℓ+12𝑝^Ω+1cos𝐼sin2ℓ2𝑎⁵𝜇³8𝑎𝑝^𝑎+1𝑝_ℎsinℓ−𝑝_𝑘cosℓ−123𝑝^ℎ𝑝_𝑘+si𝑝n^Ω𝐼𝑝^𝐼cos2ℓ43𝑝^2ℎ−𝑝^2𝑘+𝑝^2𝐼−𝑝^2Ωsin²𝐼sin2ℓ.(5.3) Double prime is omitted to simplify the notation.

For transfers between close orbits, 𝐹₁ and 𝑆₁ can be linearized around a suitable reference orbit, and an approximate first-order analytical solution can be determined. This solution is given by Δ𝑥=𝐴𝑝₀+𝐵,(5.4) where Δ𝑥=[Δ𝛼ΔℎΔ𝑘Δ𝐼ΔΩ^]^𝑇,𝑝₀ is the 5×1 vector of initial value of the adjoint variables, 𝐴 is a 5×5 symmetric matrix concerning the optimal thrust acceleration and 𝐵 is a 5×1 vector containing the perturbative effects of the oblateness of the central body. The variables 𝛼=𝑎/𝑎 and 𝑝_𝛼=𝑎𝑝_𝑎 are introduced to make the adjoint and the state vectors dimensionless, and the variables Ω^=Ωsin𝐼 and 𝑝_Ω^′=𝑝_Ω/sin𝐼 are introduced to simplify the matrix 𝐴. The adjoint vector is constant. The matrix 𝐴 and the vector 𝐵 are given by 𝐴=⎡⎢⎢⎢⎢⎢⎣𝑎^𝛼𝛼𝑎𝑎_𝛼ℎ𝑎_𝛼𝑘00_ℎ𝛼𝑎𝑎_ℎℎ𝑎_ℎ𝑘00_𝑘𝛼0𝑎0_𝑘ℎ𝑎0_𝑘𝑘𝑎0_𝐼𝐼0𝑎_𝐼Ω^′000𝑎_Ω^′_𝐼𝑎_Ω^′_Ω^′𝐵=⎡⎢⎢⎢⎣𝑏_𝛼⎤⎥⎥⎥⎥⎥⎦,(5.5)𝑏_ℎ𝑏_𝑘𝑏_𝐼𝑏_Ω^′⎤⎥⎥⎥⎦,(5.6) where 𝑎_𝛼𝛼=4𝑎⁵𝜇³𝑎_𝛼ℎℓ|||||^ℓ𝑓_ℓ0,=𝑎_ℎ𝛼=4𝑎⁵𝜇³𝑎_𝛼𝑘sinℓ|||||^ℓ𝑓_ℓ0,=𝑎_𝑘𝛼=−4𝑎⁵𝜇³𝑎_ℎℎcosℓ|||||^ℓ𝑓_ℓ0,=𝑎⁵𝜇³52ℓ+34sin2𝑎_ℎ𝑘ℓ|||||^ℓ𝑓_ℓ0,=𝑎_𝑘ℎ=−34𝑎⁵𝜇³𝑎_𝑘𝑘cos2ℓ|||||^ℓ𝑓_ℓ0,=𝑎⁵𝜇³25ℓ−34sin2𝑎_𝐼𝐼ℓ|||||^ℓ𝑓_ℓ0,=𝑎⁵𝜇³12ℓ+14sin2𝑎_𝐼Ω′ℓ|||||^ℓ𝑓_ℓ0,=𝑎_Ω^′_𝐼=−14𝑎⁵𝜇³𝑎_Ω′Ω′cos2ℓ|||||^ℓ𝑓_ℓ0,=𝑎⁵𝜇³12ℓ−14sin2𝑏_𝛼ℓ|||||^ℓ𝑓_ℓ0,𝑏_ℎ=𝜀sin²𝐼cos2ℓ|||^ℓ𝑓_ℓ0,=𝜀cosℓ+14sin²𝐼−5cosℓ+73cos3𝑏_𝑘ℓ|||^ℓ𝑓_ℓ0,=𝜀sinℓ+14sin²𝐼−7sinℓ+73sin3𝑏_𝐼ℓ|||^ℓ𝑓_ℓ0,=𝜀4sin2𝑏_Ω′𝐼cos2ℓ|||^ℓ𝑓_ℓ0,=−𝜀2sin2𝐼ℓ−12sin2ℓ|||^ℓ𝑓_ℓ0,(5.7) with ℓ_𝑓=ℓ₀+𝑛(𝑡_𝑓−𝑡₀),𝑡₀ is the initial time, 𝑡_𝑓 is the final time, and 𝜀=(3/2)𝐽₂(𝑎_𝑒/𝑎)². The overbar denotes the orbital elements of the reference orbit.

The analytical solution defined by (5.4)–(5.7) is in agreement with the ones obtained through different approaches in [6, 10]. In [6] the optimization problem is formulated with Gauss equations in nonsingular orbital elements as state equations, and in [10] Hori method is applied after the transformation of variables. Equations (5.4)–(5.7) represent a complete first-order analytical solution for optimal low-thrust limited-power transfers between close quasicircular nonequatorial orbits in a gravity field that includes the second zonal harmonic 𝐽₂ in the gravitational potential. These equations contain arbitrary constants of integration that must be determined to satisfy the two-point boundary value problem of going from the initial orbit at the time 𝑡₀=0 to the final orbit at the specified final time 𝑡_𝑓=𝑇. Since they are linear in these constants, the boundary value problem can be solved by simple techniques.

An approximate expression for the optimal thrust acceleration 𝜸^∗ can be obtained from (3.5) and (5.2), by using the expansions of the elliptic motion, and it is given, up to the zeroth order in eccentricity, by 𝜸^∗=1𝑛𝑎𝑝^ℎsinℓ−𝑝_𝑘cosℓ𝐞_𝑟+2𝑝_𝛼+𝑝_ℎcosℓ+𝑝_𝑘sinℓ𝐞_𝑠+𝑝_𝐼cosℓ+𝑝_Ω^′sinℓ𝐞_𝑤.(5.8)

According to Section 2, the optimal fuel consumption variable 𝐽 is determined by simple quadrature of the last differential equation in (2.3), with 𝜸^∗ given by (5.8). An approximate expression for 𝐽 is given by Δ𝐽=12𝑎^𝛼𝛼𝑝^2𝛼+2𝑎_𝛼ℎ𝑝_𝛼𝑝_ℎ+2𝑎_𝛼𝑘𝑝_𝛼𝑝_𝑘+𝑎_ℎℎ𝑝^2ℎ+2𝑎_ℎ𝑘𝑝_ℎ𝑝_𝑘+𝑎_𝑘𝑘𝑝^2𝑘+𝑎_𝐼𝐼𝑝^2𝐼+2𝑎_ℎΩ^′𝑝_𝐼𝑝_Ω^′+𝑎_Ω^′_Ω^′𝑝^2Ω_′.(5.9)

For transfers involving a large number of revolutions, (5.4)–(5.7) can be greatly simplified by neglecting the short periodic terms in comparison to the secular ones, and, the system of algebraic can be solved analytically. Thus, 𝑝_𝛼=14𝜇³𝑝_ℎ𝑎⁵ΔΔℓ𝛼,=25𝜇³𝑝_𝑘𝑎⁵ΔΔℓℎ,=25𝜇³𝑝_𝐼𝑎⁵ΔΔℓ𝑘,=2𝜇³𝑝_Ω′𝑎⁵ΔΔℓ𝐼,=2𝜇³𝑎⁵ΔΩsin𝐼Δℓ+2𝜀sin2𝐼.(5.10) Introducing these equations into (5.8) and (5.9), 𝜸^∗ and 𝐽 are obtained explicitly in terms of the imposed variations of the orbital elements, the perturbative effects due to the oblateness of the central body and the duration of the maneuver: 𝜸^∗=𝜇𝑎²25ΔΔℓℎsinℓ−Δ𝑘Δℓcosℓ𝐞_𝑟+12ΔΔℓ𝛼+45ΔΔℓℎcosℓ+Δ𝑘Δℓsinℓ𝐞_𝑠+2Δ𝐼Δℓcosℓ+ΔΩΔℓ+𝜀cos𝐼sin𝐼sinℓ𝐞_𝑤Δ𝐽,(5.11)Δℓ=41𝜇³𝑎⁵12ΔΔℓ𝛼²+45Δℎ²+Δ𝑘²Δℓ²+4ΔΔℓ𝐼²+ΔΩΔℓ+𝜀cos𝐼²sin²𝐼.(5.12)

5.2. Transfers between Close Quasicircular Near-Equatorial Orbits

Consider the Mathieu transformation [33] defined by 𝑃_𝑎^=𝑎^,𝜉^=𝑒^cos𝜔^+Ω^,𝜂^=𝑒^sin𝜔^+Ω^,=sin𝐼^2cosΩ^,𝑄^=sin𝐼^2sinΩ^𝑝,𝜆^=𝑀^+𝜔^+Ω^,(5.13)_𝑎𝑝_𝜉=𝑝_𝑎,=𝑝^𝑒cos𝜔^+Ω^−𝑝^𝜔−𝑝^𝑀𝑒^sin𝜔^𝑝_𝜂+Ω^,=𝑝^𝑒sin𝜔^+Ω^+𝑝^𝜔−𝑝^𝑀𝑒^cos𝜔^𝑝_𝑃+Ω^,=𝑝^𝐼2cosΩ^cos𝐼^/2−𝑝^Ω−𝑝^𝜔sinΩ^𝑝sin𝐼^/2_𝑄,=𝑝^𝐼2sinΩ^cos𝐼^/2+𝑝^Ω−𝑝^𝜔cosΩ^𝑝sin𝐼^/2_𝜆,=𝑝^𝑀.(5.14) The new set of canonical variables, designated by the double prime, is nonsingular for circular equatorial orbits.

Introducing (5.14) into (4.6) and using the expansions of the elliptic motion in terms of eccentricity and mean anomaly [31], one finds that the new averaged maximum Hamiltonian and the generating function are written up to the zeroth order in eccentricity and first-order in inclination as 𝐹₁=32𝑛𝐽²𝑎𝑎^𝑒²𝑄𝑝_𝑃−𝑃𝑝_𝑄+𝑎2𝜇4𝑎²𝑝^2𝑎+52𝑝^2𝜉+𝑝^2𝜂+18𝑝^2𝑃𝑆+𝑝^2𝑄,(5.15)₁=32𝐽²𝑎𝑎^𝑒²𝑝^𝜉cos𝜆+𝑝_𝜂sin𝜆+12𝑄𝑝^𝑃+𝑃𝑝_𝑄sin2𝜆+12𝑃𝑝^𝑃+1−𝑄𝑝_𝑄cos2𝜆2𝑎⁵𝜇³8𝑎𝑝^𝑎−3𝑝_𝜉sin𝜆−𝑝_𝜂cos𝜆+34𝑝^2𝜉−𝑝^2𝜂+116𝑝^2𝑃−𝑝^2𝑄sin2𝜆2𝑝^𝜉𝑝_𝜂+18𝑝^𝑃𝑝_𝑄cos2𝜆.(5.16) Double prime is omitted to simplify the notation.

For transfers between close orbits, 𝐹₁ and 𝑆₁ can be linearized around a suitable reference orbit, and an approximate first-order analytical solution can be determined. This solution is given by Δ𝑥=𝐶𝑝₀+𝐷,(5.17) where Δ𝑥=[Δ𝛼Δ𝜉Δ𝜂Δ𝑃Δ𝑄]^𝑇,𝑝₀ is the 5×1 vector of initial value of the adjoint variables, 𝐶 is a 5×5 matrix concerning the optimal thrust acceleration, and 𝐷 is a 5×1 vector containing the perturbative effects of the oblateness of the central body. The matrix 𝐶 and the vector 𝐷 are given by 𝐶=⎡⎢⎢⎢⎢⎢⎣𝑐^𝛼𝛼𝑐𝑐_𝛼𝜉𝑐_𝛼𝜂00_𝜉𝛼𝑐𝑐_𝜉𝜉𝑐_𝜉𝜂00_𝜂𝛼0𝑐0_𝜂𝜉𝑐0_𝜂𝜂𝑐0_𝑃𝑃0𝑐_𝑃𝑄000𝑐_𝑄𝑃𝑐_𝑄𝑄𝐷=⎡⎢⎢⎢⎣𝑑0_𝜉⎤⎥⎥⎥⎥⎥⎦,(5.18)𝑑_𝜂𝑑_𝑃𝑑_𝑄⎤⎥⎥⎥⎦,(5.19) where 𝑐_𝛼𝛼=4𝑎⁵𝜇³𝜆|||||^𝜆𝑓_𝜆0,𝑐_𝛼𝜉=𝑐_𝜉𝛼=4𝑎⁵𝜇³𝑐_𝛼𝜂sin𝜆|||||^𝜆𝑓_𝜆0,=𝑐_𝜂𝛼=−4𝑎⁵𝜇³cos𝜆|||||^𝜆𝑓_𝜆0,𝑐_𝜉𝜉=𝑎⁵𝜇³25𝜆+34sin2𝑐_𝜉𝜂𝜆|||||^𝜆𝑓_𝜆0,=𝑐_𝜂𝜉=−34𝑎⁵𝜇³cos2𝜆|||||^𝜆𝑓_𝜆0,𝑐_𝜂𝜂=𝑎⁵𝜇³25𝜆−34sin2𝑐_𝑃𝑃𝜆|||||^𝜆𝑓_𝜆0,=18𝑎⁵𝜇³Δ𝜆cos𝜀Δ𝜆+12sin(2+𝜀)𝜆_𝑓−𝜀𝜆₀−12sin(2+𝜀)𝑐𝜆₀−𝜀𝜆_𝑓_𝑃𝑄,=18𝑎⁵𝜇³Δ𝜆sin𝜀Δ𝜆−12cos(2+𝜀)𝜆_𝑓−𝜀𝜆₀+12cos(2+𝜀)𝑐𝜆₀−𝜀𝜆_𝑓_𝑄𝑃,=18𝑎⁵𝜇³−Δ𝜆sin𝜀Δ𝜆−12cos(2+𝜀)𝜆_𝑓−𝜀𝜆₀+12cos(2+𝜀)𝑐𝜆₀−𝜀𝜆_𝑓_𝑄𝑄,=18𝑎⁵𝜇³Δ𝜆cos𝜀Δ𝜆−12sin(2+𝜀)𝜆_𝑓−𝜀𝜆₀+12sin(2+𝜀)𝑑𝜆₀−𝜀𝜆_𝑓_𝜉,=32𝐽²𝑎𝑎^𝑒²cos𝜆||||^𝜆𝑓_𝜆0,𝑑_𝜂=32𝐽²𝑎𝑎^𝑒²𝑑_𝑃sin𝜆||||^𝜆𝑓_𝜆0,=𝑃₀cos𝜀Δ𝜆−12𝜀cos𝜀𝜆_𝑓−(2+𝜀)𝜆₀+12𝜀cos𝜀+𝑄₀𝜆₀−(2+𝜀)𝜆_𝑓−1𝑑sin𝜀Δ𝜆+12𝜀sin𝜀𝜆_𝑓−(2+𝜀)𝜆₀−12𝜀sin𝜀𝜆₀−(2+𝜀)𝜆_𝑓_𝑄,=𝑄₀cos𝜀Δ𝜆+12𝜀cos𝜀𝜆_𝑓−(2+𝜀)𝜆₀−12𝜀cos𝜀+𝑃₀𝜆₀−(2+𝜀)𝜆_𝑓−1−sin𝜀Δ𝜆+12𝜀sin𝜀𝜆_𝑓−(2+𝜀)𝜆₀−12𝜀sin𝜀𝜆₀−(2+𝜀)𝜆_𝑓.(5.20) The adjoint variables are 𝑝𝑝_𝑎𝜉=𝑝_𝑎0,𝑝_𝜂=𝑝_𝜉0,𝑝_𝑃=𝑝_𝜂0,𝑝=𝑝_𝑃0cos𝜀Δ𝜆+𝑝_𝑄0sin𝜀Δ𝜆,_𝑄=−𝑝_𝑃0sin𝜀Δ𝜆+𝑝_𝑄0cos𝜀Δ𝜆,(5.21) where Δ𝜆=𝜆_𝑓−𝜆₀=𝑛(𝑡_𝑓−𝑡₀). The overbar denotes the orbital elements of the reference orbit. Equations (5.20) do not contain higher order terms in 𝜀  as appear implicitly in the solution presented in [7].

The analytical solution defined by (5.17)–(5.21) is in agreement with the ones obtained through different approaches in [7, 10]. In [7] the optimization problem is formulated with Gauss equations in nonsingular orbital elements as state equations, and in [10] Hori method is applied after the transformation of variables. Equations (5.17)–(5.21) represent a complete first-order analytical solution for optimal low-thrust limited-power transfers between close quasicircular near-equatorial orbits in a gravity field that includes the second zonal harmonic 𝐽₂ in the gravitational potential. These equations contain arbitrary constants of integration that must be determined to satisfy the two-point boundary value problem of going from the initial orbit at the time 𝑡₀=0 to the final orbit at the specified final time 𝑡_𝑓=𝑇. Since they are linear in these constants, the boundary value problem can be solved by simple techniques.

As described in Section 5.1, approximate expressions can be obtained for the optimal thrust acceleration 𝜸^∗ and the fuel consumption 𝐽. These expressions are given by 𝜸^∗=1𝑛𝑎𝑝^𝜉sin𝜆−𝑝_𝜂cos𝜆𝐞_𝑟+2𝑝_𝛼+𝑝_𝜉cos𝜆+𝑝_𝜂sin𝜆𝐞_𝑠+12𝑝^𝑃Δ𝐽=1cos𝜆+𝑝_𝑄sin𝜆𝐞_𝑤,2𝑐^𝛼𝛼𝑝^2𝛼+2𝑐_𝛼𝜉𝑝_𝛼𝑝_𝜉+2𝑐_𝛼𝜂𝑝_𝛼𝑝_𝜂+𝑐_𝜉𝜉𝑝^2𝜉+2𝑐_𝜉𝜂𝑝_𝜉𝑝_𝜂+𝑐_𝜂𝜂𝑝^2𝜂+𝑐_𝑃𝑃𝑝^2𝑃+𝑐_𝑃𝑄+𝑐_𝑄𝑃𝑝_𝑃𝑝_𝑄+𝑐_𝑄𝑄𝑝^2𝑄.(5.22)

For transfers involving a large number of revolutions, (5.17)–(5.21) can be greatly simplified by neglecting the short periodic terms in comparison to the secular ones, and, the system of algebraic can be solved analytically. Thus, 𝑝_𝛼=14𝜇³𝑝_𝜉𝑎⁵ΔΔ𝛼𝜆,=25𝜇³𝑝_𝜂𝑎⁵ΔΔ𝜆𝜉,=25𝜇³𝑝_𝑃𝑎⁵ΔΔ𝜆𝜂,=𝜇³𝑎⁵8Δ𝜆𝑃⁰𝑝_𝑄+Δ𝑃cos𝜀Δ𝜆−𝑄₀+Δ𝑄sin𝜀Δ𝜆−𝑃₀,=𝜇³𝑎⁵8Δ𝜆𝑄⁰+Δ𝑄cos𝜀Δ𝜆+𝑃₀+Δ𝑃sin𝜀Δ𝜆−𝑄₀.(5.23) The subscript denoting the constants is omitted. Introducing these equations into (5.22), 𝜸^∗ and 𝐽 are obtained explicitly in terms of the imposed variations of the orbital elements, the perturbative effects due to the oblateness of the central body and the duration of the maneuver: 𝜸^∗=𝜇𝑎²25ΔΔ𝜆𝜉sin𝜆−Δ𝜂Δ𝜆cos𝜆𝐞_𝑟+12ΔΔ𝛼𝜆+45ΔΔ𝜆𝜉cos𝜆+Δ𝜂Δ𝜆sin𝜆𝐞_𝑠+4Δ𝜆𝑃^𝑓cos𝜆+𝜃_𝑓−𝑃⁰cos𝜆+𝜃₀+𝑄^𝑓sin𝜆+𝜃_𝑓−𝑄⁰sin𝜆+𝜃₀𝐞^𝑤Δ𝐽,Δ𝜆=14𝜇³𝑎⁵12ΔΔ𝛼𝜆²+45Δ𝜉²+Δ𝜂²+16cos𝜀ΔΔ𝜆²𝜆Δ𝑃²+Δ𝑄²Δ𝜆²+Δ4𝜆²sin𝜀Δ𝜆×𝑃_𝑓2𝑃⁰sin𝜀Δ𝜆2−𝑄⁰cos𝜀Δ𝜆+𝑄_𝑓2𝑃⁰cos𝜀Δ𝜆2+𝑄⁰sin𝜀Δ𝜆2,(5.24) with 𝜃₀=𝜀𝑛(𝑡−𝑡₀) and 𝜃_𝑓=𝜀𝑛(𝑡−𝑡_𝑓).

5.3. General Remarks on Optimal Transfers between Close Quasicircular Orbits

Considering the analytical solutions described above, general remarks about the effects of the oblateness of the central body on optimal transfers between close quasicircular orbits can be stated.

(1) Matrices 𝐴 and 𝐶, defined by (5.5) and (5.18), respectively, can be decomposed into two square matrices 3×3 and 2×2. This fact shows that there exists an uncoupling between the in-plane modifications and the rotation of the orbital plane. Similar results were obtained by Edelbaum [1] and Marec [3] for transfers in a Newtonian central field.(2) There is a normal component of the optimal thrust acceleration that counteracts the perturbative effects due to the oblateness of the central body. This normal component is proportional to 𝜀.(3) In general, the maneuvers in a noncentral gravity field are more expensive than the maneuvers in Newtonian central field, taking into account that the perturbative effects due to the oblateness of the central body must be counteracted.(4) Fuel can be saved for long-time maneuvers which involves modification of the longitude of the ascending node if the terminal orbits are direct (0^∘<𝐼<90^∘) and ΔΩ<0, and if the terminal orbits are retrograde (90^∘<𝐼<180^∘) and ΔΩ>0.(5) The extra consumption needed to counteract the perturbative effects due to the oblateness of the central body reaches its maximum for maneuvers between orbits with 45^∘ or 135^∘ of inclination and its minimum (null extra consumption) for maneuvers between equatorial or polar orbits (the gravity field is symmetric in these cases). (6) As the transfer time increases, the extra fuel consumption increases.(7) The extra fuel consumption is greater for transfers between low orbits.

6. Optimal Long-Time Transfers between Arbitrary Orbits

In this section the two-point boundary-value problem for long-time transfers between arbitrary orbits is formulated. In Section 7, this boundary-problem problem will be solved numerically through a shooting method, and the solution will be compared to the analytical ones in the case of orbits with small eccentricities.

Since the averaged maximum Hamiltonian, given by (4.6), becomes singular for circular and/or equatorial orbits, a set of nonsingular elements is introduced.

Consider the Mathieu transformation [33] defined by 𝜉_𝑎^=𝑎^,𝜂_=𝑒^cos𝜔^+Ω^,𝑃_=𝑒^sin𝜔^+Ω^,=sin𝐼^2cosΩ^𝑄_,=sin𝐼^2sinΩ^𝑝,(6.1)_𝑎𝑝_𝜉=𝑝^𝑎,=𝑝^𝑒cos𝜔^+Ω^−𝑝^𝜔−𝑝^𝑀𝑒^sin𝜔^𝑝_𝜂+Ω^,=𝑝^𝑒sin𝜔^+Ω^+𝑝^𝜔−𝑝^𝑀𝑒^cos𝜔^𝑝_𝑃+Ω^,=𝑝^𝐼2cosΩ^cos𝐼^/2−𝑝^Ω−𝑝^𝜔sinΩ^𝑝sin𝐼^/2_𝑄,=𝑝^𝐼2sinΩ^cos𝐼^/2+𝑝^Ω−𝑝^𝜔cosΩ^sin𝐼^/2.(6.2) This transformation of variables is the same one given by (5.12), excepting the fast-phase, which is unnecessary (short-periodic terms have been eliminated). Double prime designates the new variables.

Introducing (6.2) into (4.6), 𝐹₁ is written as 𝐹₁=−32𝑛𝐽²𝑎𝑎^𝑒²1−𝜉²−𝜂²+1⁻²1−2𝑃²+𝑄²𝜉𝑝_𝜂−𝜂𝑝_𝜉+𝑃𝑝_𝑄−𝑄𝑝_𝑃21−51−2𝑃²+𝑎+𝑄²²𝜉𝑝^𝜂−𝜂𝑝_𝜉2𝜇4𝑎²𝑝^2𝑎+𝑝^2𝜉52−25𝜉²−2𝜂²+𝑝^2𝜂52−52𝜂²−2𝜉²−𝜉𝜂𝑝^𝜉𝑝_𝜂+121−𝜉²−𝜂²⁻¹1+32𝜉²+3×12𝜂²4𝑝^2𝑃+𝑝^2𝑄−14𝑃𝑝^𝑃+𝑄𝑝_𝑄²+141−𝑃²−𝑄²⁻¹×+4𝑄𝑝𝑃_𝑃2−𝑃𝑝_𝑄²+4𝜉𝑝_𝜂−𝜂𝑝_𝜉𝑃𝑝_𝑄−𝑄𝑝_𝑃−5+𝑄²𝜉𝑝_𝜂−𝜂𝑝_𝜉²81−𝑃²×2𝜉𝑝_𝜂−𝑄²⁻¹𝜉²−𝜂²𝑃²−𝑄²+4𝜉𝜂𝑃𝑄−5−𝜂𝑝_𝜉+𝑃𝑝_𝑄−𝑄𝑝_𝑃²+1−𝑃²−𝑄²𝑝^2𝑃+𝑝^2𝑄2𝜉𝜂𝑝^𝑃𝑝_𝑄+58𝜉²+5−𝜂²𝑝^2𝑃−𝑝^2𝑄2𝜉𝑝^𝜂−𝜂𝑝_𝜉−𝜉²−𝜂²𝑃𝑝_𝑄+𝑄𝑝_𝑃+2𝜉𝜂𝑃𝑝_𝑃−𝑄𝑝_𝑄.(6.3) Double prime is omitted to simplify the notation. Note that (5.15) is a simplification of (6.3) for quasicircular near-equatorial orbits.

The two-point boundary-value problem for long-time transfers between arbitrary orbits is formulated by the following system of canonical differential equations: 𝑑𝑥𝑑𝑡=𝜕𝜕𝐹𝑝¹^𝑇,𝑑𝑝𝑑𝑡=−𝜕𝜕𝐹𝑥¹^𝑇,(6.4) where 𝑥=[𝑎𝜉𝜂𝑃𝑄]^𝑇 and 𝑝=[𝑝_𝑎𝑝_𝜉𝑝_𝜂𝑝_𝑃𝑝_𝑄]^𝑇, and the boundary-conditions 𝑎𝑎𝑡_𝑓𝑡₀=𝑎₀,𝜉𝑡₀=𝜉₀,𝜂𝑡₀=𝜂₀,𝑃𝑡₀=𝑃₀,𝑄𝑡₀=𝑄₀,=𝑎_𝑓,𝜉𝑡_𝑓=𝜉_𝑓,𝜂𝑡_𝑓=𝜂_𝑓,𝑃𝑡_𝑓=𝑃_𝑓,𝑄𝑡_𝑓=𝑄_𝑓,(6.5) with initial time 𝑡₀=0 and final time 𝑡_𝑓=𝑇. The explicit form of state and adjoint equations is derived by using MAPLE Software 9, and they are not presented in text because they are extensive.

7. Results

In this section, numerical and analytical results are compared for maneuvers described in Tables 1 and 2. Earth is the central body, such that 𝐽₂=1.0826×10⁻³. In Table 1, orbital elements of initial and final orbits are defined for a maneuver of modification of all orbital elements considering different values of transfer durations 𝑇=𝑡_𝑓−𝑡₀ and ratios 𝑎_𝑒/𝑎₀. In Table 2, orbital elements of initial orbit, values of transfer duration and ratio 𝜌=𝑦_𝑓/𝑦₀, where 𝑦 denotes the orbital element that is changed in the maneuver, are defined for small (orbit corrections) and moderate amplitude maneuvers of modification of semimajor axis and inclination of the orbital plane. Results are presented in Figures 3 through 8 using canonical units. In these units, 𝑎₀=1.0 and 𝜇=1.0. For example, in metric units, taking 𝑎_𝑒=6378.2km,𝜌=0.950 corresponds to 𝑎₀=6713.89km, and 𝑇=100 (time units) corresponds to 24.20 hours (1.008 day), and 𝜌=0.150 corresponds to 𝑎₀=42521.34km, and 𝑇=250 (time units) corresponds to 964.45 hours (40.185 days). Thus, maneuvers described in the tables below have duration of 1 to 74 days, approximately. On the other hand, considering that the value of |𝜸^∗| is approximately given by √2Δ𝐽/𝑇 (short periodic terms are neglected), |𝜸^∗| has values between 1.0×10⁻⁴ and 4.0×10⁻³ (canonical units) for all results presented in Figures 3–5. These values of magnitude of the thrust acceleration are typical for LP systems [5, 27].

Table 1 

Maneuver of modification of all orbital elements.

Table 2 

Small and moderate amplitude maneuvers.

Figure 3 

Consumption for maneuver of modification of all orbital elements.

(a)

(b)

(c)

(d)

Figure 4 

Consumption for maneuver of modification of semimajor axis.

(a)

(b)

(c)

(d)

Figure 5 

Consumption for maneuver of modification of inclination.

Figure 3 shows the consumption 𝐽 as function of ratio 𝑎_𝑒/𝑎₀ and values of 𝑇 for the maneuver of modification of all orbital elements described in Table 1. Figures 4 and 5 show the consumption 𝐽 as function of ratio 𝜌=𝑦_𝑓/𝑦₀ and 𝑇 for the small and moderate amplitude maneuvers described in Table 2, considering that the maneuvers are performed in central or noncentral gravity field. In these figures, solid line represents analytical solution, and dashed line represents numerical solution. Analytical solutions, including short periodic terms, are obtained by solving the linear system of algebraic equations given by (5.4)–(5.7) (transfer between nonequatorial orbits) with the semimajor axis and the inclination of the reference orbit defined by 𝑎=(𝑎₀+𝑎_𝑓)/2 and 𝐼=(𝐼₀+𝐼_𝑓)/2, respectively. In order to include the short periodic terms, it is assumed that the initial position of the vehicle is the pericenter of the initial orbit. Variables 𝐼 and Ω are transformed into nonsingular variables 𝑃 and 𝑄 to compare the analytical results with the numerical results. Numerical solutions are obtained by solving the two-point boundary value problem described by (6.4) and (6.5) through a shooting method [34].

Considering the maneuvers described in Tables 1 and 2, Figures 3–5 show that analytical and numerical solutions yield quite similar values of the consumption 𝐽 in the following cases: (i) maneuvers of modification of all orbital elements with ratio 𝑎_𝑒/𝑎₀<0.550; (ii) maneuvers of modification of semimajor axis with ratios 𝜌<1.25 and 𝑎_𝑒/𝑎₀<0.550; (iii) maneuvers of modification of inclination with ratios 𝜌<1.50 and 𝑎_𝑒/𝑎₀<0.550, for all values of 𝑇. For maneuvers with very high ratio 𝑎_𝑒/𝑎₀, analytical and numerical solutions only yield similar results for the same maneuvers described above for 𝑇=100 (time units).

Figures 6–8 represent the time history of state variables (nonsingular orbital elements) for maneuvers of modification of all orbital elements (Table 1) considering three values of ratios 𝑎_𝑒/𝑎₀—0.150, 0.550, and 0.950—and two values of transfer duration 𝑇—100 and 250 (time units). According to these figures, the analytical solution yields a good representation of the numerical solution for all state variables, excepting the semimajor axis for the maneuver with 𝑎_𝑒/𝑎₀=0.950 and 𝑇=250. The numerical solution shows that the semimajor axis is strongly perturbed by the oblateness of earth, and its time evolution cannot be described by the linear approximation given by the analytical solution. On the other hand, the analytical solutions shows that the amplitude of the short periodic terms decreases as the transfer duration increases such that these terms can be omitted for very long-time transfers.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 6 

Time history of state variables for maneuver of modification of all orbital elements with 𝑎𝑒/𝑎0=0.950---𝑇=100 and 250.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 7 

Time history of state variables for maneuver of modification of all orbital elements with 𝑎𝑒/𝑎0=0.550---𝑇=100 and 250.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 8 

Time history of state variables for maneuver of modification of all orbital elements with 𝑎𝑒/𝑎0=0.150---𝑇=100 and 250.

Figure 9 

Vectors 𝐞,𝐡, and 𝐍.

8. Conclusion

In this paper an approach based on canonical transformations is presented for the problem of optimal low-thrust limited-power maneuvers between close orbits with small eccentricities in noncentral gravity which includes the perturbative effects of the second zonal harmonic of the gravitational potential. Analytical first-order solutions are derived through this approach for transfers between nonequatorial orbits and for transfers between near-equatorial orbits. These analytical solutions have been compared to a numerical solution for long-time transfers between arbitrary orbits described by averaged maximum Hamiltonian, expressed in nonsingular orbital elements. The results show that the analytical solutions yield good estimates of the fuel consumption for preliminary mission analysis considering maneuvers with moderate amplitude, ratio 𝑎_𝑒/𝑎₀ and transfers duration. Numerical results also show that the semimajor axis is strongly perturbed by the oblateness of earth, and its time evolution cannot be described by the linear approximation given by the analytical solution.

Appendix

The general solution of the differential equations for the adjoint variables 𝐩_𝑟 and 𝐩_𝑣 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 below.

Let us consider the inverse of the point transformation defined by (3.2) and (3.3): 𝑎=𝑟𝑒2−𝑟𝑣²/𝜇₂,=1−ℎ²cos𝐼=𝐤𝜇•𝑎𝐡,cosΩ=𝐢•ℎ𝐍,cos𝜔=𝐍𝑁•𝐞,cos𝐸=1𝑁𝑒,𝑒1−𝑎𝑟,(A.1) where the eccentricity vector 𝐞, the angular momentum vector 𝐡, and the nodal vector 𝐍, shown in Figure 9, are given, as function of the Cartesian elements 𝐫 and 𝐯, by the following equations: 𝐞=1𝜇𝑣²−𝜇𝐍𝐡==𝑟𝐤𝐫××𝐫𝐯𝐡−,.(𝐫•𝐯)𝐯,(((AAA...234))) Here the symbol × denotes the cross product. Note that the true anomaly 𝑓 has been replaced by eccentric anomaly 𝐸. These anomalies are related through the following equation: tan𝑓2=1+𝑒1−𝑒tan𝐸2.(A.5)

Now, consider the variations in the Cartesian elements, 𝛿𝐫 and 𝛿𝐯, given in the moving frame of reference by 𝛿𝛿𝐯𝐫==𝛿𝛿𝑢𝜉𝐞𝐞_𝐫𝐫+𝛿𝜂𝐞_𝑠+𝛿𝜁𝐞_𝑤,+𝛿𝑣𝐞_𝑠+𝛿𝑤𝐞_𝑤.(A.6) The variations of the orbital elements—𝑎,𝑒,𝐼,Ω,𝜔, and 𝐸—induced by the variations in the Cartesian elements, 𝛿𝐫 and 𝛿𝐯, are obtained straightforwardly from (A.1) through (A.6) and are given by 𝛿𝑎=2𝑎𝑟²𝐫•𝛿𝐫𝛿𝑒=−𝐡•𝛿𝑟𝐡+2𝜇𝑎²𝐯•𝛿𝐯,𝛿𝐼=𝐡•𝛿𝐡𝜇𝑎𝑒+2𝜇ℎ𝑎²²𝑒𝛿𝑎,ℎ²cot𝐼−𝐤•𝛿𝐡𝛿Ω=𝐍•𝛿𝐍ℎcsc𝐼,𝑁²cotΩ−𝐢•𝛿𝐍𝛿𝜔=𝛿𝑒𝑁cscΩ𝛿𝑒𝐸+=𝐍𝑁•1𝛿²𝐍cot𝜔−𝑁1𝑒[𝐍•𝛿𝐞+𝐞•𝛿𝐍]csc𝜔,𝑒sin𝐸𝐫•𝑟𝑎𝛿𝐫−𝑎𝑟²𝛿𝑎+cos𝐸𝛿𝑒,(A.7) where the variations of the vectors 𝐞,𝐡, and 𝐍 are written as 𝜇𝛿𝐞=𝑣^2𝑠𝛿𝜉−𝑣_𝑠𝑣_𝑟𝛿𝜂+2𝑟𝑣_𝑠𝛿𝑣𝐞_𝐫+−𝑣_𝑠𝑣_𝑟𝛿𝜉+𝑣^2𝑟−𝜇𝑟𝛿𝜂−𝑟𝑣^𝑠𝛿𝑢−𝑟𝑣_𝑟𝛿𝑣𝐞_𝑠+𝑣²−𝜇𝑟𝛿𝜁−𝑟𝑣^𝑟𝛿𝐡=−𝑣_𝑠𝛿𝑤𝐞_𝑤,𝛿𝐍=𝑣_𝑠𝛿𝜁𝐞_𝑟+𝑣_𝑟𝛿𝜁−𝑟𝛿𝑤𝐞_𝑠+𝑣_𝑠𝛿𝜉−𝑣_𝑟𝛿𝜂+𝑟𝛿𝑣𝐞_𝐰,𝛿𝜉−𝑣_𝑟𝛿𝜂+𝑟𝛿𝑣sin𝐼cos(𝜔+𝑓)−𝑣_𝑟𝛿𝜁−𝑟𝛿𝑤cos𝐼𝐞_𝐫+−𝑣_𝑠𝛿𝜉−𝑣_𝑟𝛿𝜂+𝑟𝛿𝑣sin𝐼sin(𝜔+𝑓)−𝑣_𝑠cos𝐼𝛿𝜁𝐞_𝑠+𝑣_𝑟𝛿𝜁−𝑟𝛿𝑤sin𝐼sin(𝜔+𝑓)+𝑣_𝑠sin𝐼cos(𝜔+𝑓)𝛿𝜁𝐞_𝑤.(A.8) Here 𝑣_𝑟 and 𝑣_𝑠 denote the radial and circumferential components of the velocity vector, respectively (see (3.3).

In the moving fame of reference, 𝐞,𝐡,𝐍, and the unit vectors 𝐢 and 𝐤 are written as 𝐞=ℎ𝑣^𝑠𝜇−1𝐞^𝐫−ℎ𝑣^𝑟𝜇𝐞^𝑠𝐡=ℎ𝐞_𝐰,(A.9)𝐍=ℎsin=𝐼cos𝜇𝑎(𝜔1+−𝑓𝑒)²𝐞_𝐫𝐞_𝐰,(A.10)𝐢=(cosΩcos(𝜔+𝑓)−s−inℎΩsisnin𝐼(s𝜔in+(𝜔𝑓)+co𝑓s)𝐼𝐞)_𝑠,𝐞_𝑟(A.11)𝐤=−sin(co(𝜔sΩ+s𝑓in)s(i𝜔n𝐼+𝐞𝑓_𝑟)+sinΩcos(𝜔+𝑓)cos𝐼)𝐞_𝑠+sinΩsin𝐼𝐞_𝑤,(A.12)+cos(𝜔+𝑓)sin𝐼𝐞_𝑠+cos𝐼𝐞_𝑤.(A.13)

From (3.2), (3.3), and (A.2) through (A.11), one gets the explicit form of the variations of the orbital elements—𝑎,𝑒,𝐼,Ω,𝜔, and 𝐸—induced by the variations in the Cartesian elements, 𝛿𝐫 and 𝛿𝐯. At this point, it is useful to replace the variation in the eccentric anomaly 𝐸 by the variation in the mean anomaly 𝑀, obtained from the well-known Kepler’s equation 𝑀=𝐸−𝑒sin𝐸.(A.14) This variation is given by 𝛿𝑀=𝑟𝑎𝛿𝐸−sin𝐸𝛿𝑒.(A.15) Thus, 𝛿𝑎=2𝑎𝑟²𝛿𝜉+2𝑒sin𝑓𝑛√1−𝑒²𝛿𝑢+2√1−𝑒²𝛿𝑒=𝑎1−𝑒₂𝑛𝑎𝑟𝛿𝑣,𝑟²cos𝐸𝛿𝜉+sin𝑓𝑎𝛿𝜂+√1−𝑒²𝑛𝑎sin𝑓𝛿𝑢+√1−𝑒²𝛿𝐼=1𝑛𝑎(cos𝐸+cos𝑓)𝛿𝑣,𝛿Ω=𝑟1cos𝐸sin𝜔+si√n𝐸1−co𝑒s₂𝜔𝛿𝜁+𝑛𝑎√11−𝑒²𝑎𝑟cos(𝜔+𝑓)𝛿𝑤,𝛿𝜔=𝑟sisnin𝑓𝐼−cos𝐸cos𝜔+si√n𝐸1−sin𝑒²𝜔𝛿𝜁+𝑛𝑎sin𝐼1√1−𝑒²𝑎𝑟sin(𝜔+𝑓)𝛿𝑤,−𝑒√𝑟𝛿𝜉−𝑎𝑒𝑒+1c−os𝑒𝑓²𝛿𝜂+co𝑟t𝐼cos𝐸cos𝜔−si√n𝐸1−sin𝑒²𝜔𝛿𝜁1−𝑒²𝑛𝑎𝑒cos𝑓𝛿𝑢+√1−𝑒²−cot𝐼𝑛𝑎𝑒sin𝑓1+1+𝑒1cos𝑓𝛿𝑣𝑛𝑎√1−𝑒²𝛿𝑀=−1−𝑒₃𝑎𝑟sin(𝜔+𝑓)𝛿𝑤,cos𝐸𝑒𝑟√1−𝑒²sin𝑓𝛿𝜉+√1−𝑒²+1−𝑒₂𝑎𝑒cos𝑓𝛿𝜂𝑛𝑎𝑒cos𝑓−1+𝑒2𝑒cos𝑓𝛿𝑢−1−𝑒²𝑛𝑎𝑒sin𝑓1+1+𝑒1cos𝑓𝛿𝑣.(A.16)

Equations (A.16) can be put in the form ⎡⎢⎢⎢⎢⎣𝛿𝛿𝛿𝛿𝛿𝛿𝑀Ω𝜔𝐼𝑎𝑒⎤⎥⎥⎥⎥⎦=Δ₋₁⎡⎢⎢⎢⎢⎣𝛿𝛿𝛿𝛿𝛿𝛿𝑤𝜁𝑣𝜂𝜉𝑢⎤⎥⎥⎥⎥⎦,(A.17) where the matrix Δ⁻¹ 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 [32], the general solution of the differential equations for the adjoint variables 𝐩_𝑟 and 𝐩_𝑣 is given by 𝐩𝐩^𝐫_𝐯=Δ⁻¹^𝑇⎡⎢⎢⎢⎢⎣𝑝𝑝^𝑎_𝑒𝑝_𝐼𝑝_Ω𝑝_𝜔𝑝_𝑀⎤⎥⎥⎥⎥⎦,(A.18) with 𝐩_𝑟 and 𝐩_𝑣 expressed in the moving frame of reference. Equations (3.4) and (3.5) are obtained straightforwardly from the above equation.

Acknowledgment

This research has been supported by CNPq under contract 305049/2006-2.