Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2009 / Article

Research Article | Open Access

Volume 2009 |Article ID 503168 | 35 pages | https://doi.org/10.1155/2009/503168

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

Academic Editor: Dane Quinn
Received12 Oct 2008
Revised14 Feb 2009
Accepted15 Mar 2009
Published19 Jul 2009

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 𝐽2 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 𝐽2 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/𝑔0, is between 10βˆ’4 and 10βˆ’2. For such system, the fuel consumption is described by the variable 𝐽 defined as 1𝐽=2ξ€œπ‘‘π‘‘0𝛾2𝑑𝑑,(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π‘š0ξ‚Ά,(2.2) where 𝑃max is the maximum power, and π‘š0 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.

The optimization problem is formulated as follows: it is proposed to transfer the space vehicle 𝑀 from the initial state (𝐫0,𝐯0,0) at the initial time 𝑑0=0 to the final state (𝐫𝑓,𝐯𝑓,𝐽𝑓) at the specified final time 𝑑𝑓, such that the final consumption variable 𝐽𝑓 is a minimum. The state equations are 𝑑𝐫𝑑𝑑=𝐯,𝑑𝐯𝑑𝑑=𝐠(𝐫)+𝜸,𝑑𝐽=1𝑑𝑑2𝛾2,(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 π‘‘π‘“βˆ’π‘‘0 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), 𝐻=π©π‘Ÿβ€’π―+𝐩𝑣1β€’(𝐠(𝐫)+𝜸)+2𝑝𝐽𝛾2,(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] πœ‡π (𝐫)=βˆ’π‘Ÿ3𝐫+πœ•π‘‰2ξ‚Άπœ•π«π‘‡,(2.9) with the disturbing function 𝑉2 defined by 𝑉2=βˆ’πœ‡π½2ξ‚΅π‘Žπ‘’2π‘Ÿ3𝑃2(sinπœ‘),(2.10) where πœ‡ is the gravitational parameter, π‘Žπ‘’ is the mean equatorial radius of the central body, 𝐽2 is the coefficient for the second zonal harmonic of the potential, πœ‘is the latitude, and 𝑃2 is the Legendre polynomial of order 2. (πœ•π‘‰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: π»βˆ—=𝐻0+𝐻𝐽2+π»π›Ύβˆ—,(2.11) where 𝐻0=π©π‘Ÿβ€’π―βˆ’π©π‘£β€’πœ‡π‘Ÿ3𝐻𝐫,𝐽2=π©π‘£β€’ξ‚΅πœ•π‘‰2ξ‚Άπœ•π«π‘‡,π»π›Ύβˆ—=𝑝𝑣22.(2.12)𝐻0 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 𝐻0: 𝑑𝐫𝑑𝑑=𝐯,π‘‘π―πœ‡π‘‘π‘‘=βˆ’π‘Ÿ3𝐫,𝑑𝐩𝐫=πœ‡π‘‘π‘‘π‘Ÿ3ξ€·π©π‘£ξ€·π©βˆ’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βˆ’π‘’2ξ€Έπž1+𝑒cos𝑓𝐫,ξƒŽ(3.2)𝐯=πœ‡π‘Žξ€·1βˆ’π‘’2𝑒sinπ‘“πžπ«+(1+𝑒cos𝑓)πžπ¬ξ€»,𝐩(3.3)𝐫=π‘Žπ‘Ÿ2ξƒ―2π‘Žπ‘π‘Ž+ξ€·ξ€·1βˆ’π‘’2𝑝cos𝐸𝑒+ξ‚€π‘Ÿπ‘Žξ‚sinπ‘“π‘’ξƒ©π‘πœ”βˆ’ξ€·1βˆ’π‘’3ξ€Έcos𝐸√1βˆ’π‘’2π‘π‘€πžξƒͺ𝐫+ξƒ―sinπ‘“π‘Žπ‘π‘’βˆ’(𝑒+cos𝑓)ξ€·π‘Žπ‘’1βˆ’π‘’2ξ€Έπ‘πœ”+√1βˆ’π‘’2cosπ‘“π‘π‘Žπ‘’π‘€ξƒ°πžπ¬+1π‘Žβˆš1βˆ’π‘’2π‘Žξ‚†ξ‚€π‘Ÿξ‚ξ‚ƒπ‘sin𝐸𝐼𝑝cosπœ”+Ξ©sinπΌβˆ’π‘πœ”ξ‚ξ‚„+√cot𝐼sinπœ”1βˆ’π‘’2ξ‚€π‘Žπ‘Ÿξ‚ξ‚ƒπ‘cos𝐸𝐼𝑝sinπœ”βˆ’Ξ©sinπΌβˆ’π‘πœ”ξ‚πžcot𝐼cosπœ”ξ‚„ξ‚‡π°,𝐩(3.4)𝐯=1βˆšπ‘›π‘Ž1βˆ’π‘’2ξƒ―ξƒ―2π‘Žπ‘’sinπ‘“π‘π‘Ž+ξ€·ξ€·1βˆ’π‘’2𝑝sinπ‘“π‘’βˆ’ξ€·1βˆ’π‘’2ξ€Έcosπ‘“π‘’π‘πœ”+ξ€·1βˆ’π‘’2ξ€Έ3/2𝑒cosπ‘“βˆ’2𝑒𝑝1+𝑒cosπ‘“π‘€ξƒ°πžπ«+ξƒ―ξ€·2π‘Ž1βˆ’π‘’2ξ€Έξ‚€π‘Žπ‘Ÿξ‚π‘π‘Ž+ξ€·1βˆ’π‘’2ξ€Έ(cos𝑓+cos𝐸)𝑝𝑒+ξ€·1βˆ’π‘’2ξ€Έsin𝑓𝑒11+𝑝1+𝑒cosπ‘“πœ”βˆ’βˆš1βˆ’π‘’2π‘π‘€ξ‚ξƒ°πžπ¬+π‘Ÿξ‚†ξ‚€π‘Žξ‚cos(πœ”+𝑓)𝑝𝐼+ξ‚€π‘Ÿπ‘Žξ‚ξ‚€π‘sin(πœ”+𝑓)Ξ©sinπΌβˆ’π‘πœ”πž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, 𝑛=βˆšπœ‡/π‘Ž3 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: 𝑓tan2=ξ‚™1+𝑒𝐸1βˆ’π‘’tan2,𝑀=πΈβˆ’π‘’sin𝐸.(3.6) The unit vectors 𝐞𝐫,𝐞𝐬, and 𝐞𝐰 of the moving frame of reference are written in the fixed frame of reference as 𝐞𝐫+𝐞=(cosΞ©cos(πœ”+𝑓)βˆ’sinΞ©sin(πœ”+𝑓)cos𝐼)𝐒(sinΞ©cos(πœ”+𝑓)+cosΞ©sin(πœ”+𝑓)cos𝐼)𝐣+sin(πœ”+𝑓)sin𝐼𝐀,𝐬𝐞=βˆ’(cosΞ©sin(πœ”+𝑓)+sinΞ©cos(πœ”+𝑓)cos𝐼)𝐒+(βˆ’sinΞ©sin(πœ”+𝑓)+cosΞ©cos(πœ”+𝑓)cos𝐼)𝐣+cos(πœ”+𝑓)sin𝐼𝐀,𝐰=sinΞ©sinπΌπ’βˆ’cosΞ©sin𝐼𝐣+cos𝐼𝐀.(3.7)

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 π»βˆ—=𝐻0+𝐻𝐽2+π»π›Ύβˆ—,(3.8) where 𝐻0=𝑛𝑝𝑀,𝐻𝐽2=2π‘›π‘Žπœ•π‘‰ξ…ž2π‘πœ•π‘€π‘Ž+√1βˆ’π‘’2π‘›π‘Ž2π‘’ξ‚Έβˆ’πœ•π‘‰ξ…ž2+βˆšπœ•πœ”1βˆ’π‘’2πœ•π‘‰ξ…ž2ξ‚Ήπ‘πœ•π‘€π‘’+1π‘›π‘Ž2√1βˆ’π‘’2ξ‚Έβˆ’sinπΌπœ•π‘‰ξ…ž2πœ•Ξ©+cosπΌπœ•π‘‰ξ…ž2ξ‚Ήπ‘πœ•πœ”πΌ+1π‘›π‘Ž2√1βˆ’π‘’2sinπΌπœ•π‘‰ξ…ž2π‘πœ•πΌΞ©+√1βˆ’π‘’2π‘›π‘Ž2π‘’ξƒ¬πœ•π‘‰ξ…ž2βˆ’πœ•π‘’π‘’cot𝐼1βˆ’π‘’2ξ€Έπœ•π‘‰ξ…ž2ξƒ­π‘πœ•πΌπœ”+1ξƒ¬π‘›π‘Žβˆ’2πœ•π‘‰ξ…ž2βˆ’ξ€·πœ•π‘Ž1βˆ’π‘’2ξ€Έπ‘Žπ‘’πœ•π‘‰ξ…ž2ξƒ­π‘πœ•π‘’Mξƒ°,π»βˆ—π›Ύ=12𝑛2π‘Ž2ξ€·1βˆ’π‘’212ξ€Ί(1βˆ’cos2𝑓)2π‘Žπ‘’π‘π‘Ž+ξ€·1βˆ’π‘’2𝑝𝑒2ξ€·+21βˆ’π‘’2sin2π‘“βˆ’π‘Žπ‘π‘Žπ‘πœ”βˆ’ξ€·1βˆ’π‘’2𝑝2π‘’π‘’π‘πœ”ξƒ­ξ€·+41βˆ’π‘’2ξ€Έ3/2ξ‚΅sinπ‘“βˆ’2𝑒1+𝑒cos𝑓+cosπ‘“π‘Žπ‘π‘Žπ‘π‘€+ξ€·1βˆ’π‘’2𝑝2𝑒𝑒𝑝𝑀+ξ€·1βˆ’π‘’2ξ€Έ22𝑒2(1+cos2𝑓)π‘πœ”2βˆ’2ξ€·1βˆ’π‘’2ξ€Έ5/2𝑒2ξ‚΅βˆ’2𝑒1+𝑒cos𝑓+cos𝑓cosπ‘“π‘πœ”π‘π‘€+ξ€·1βˆ’π‘’2ξ€Έ3𝑒2ξ‚΅βˆ’2𝑒1+𝑒cos𝑓+cos𝑓2𝑝𝑀2+4π‘Ž2ξ€·1βˆ’π‘’2ξ€Έ2ξ‚€π‘Žπ‘Ÿξ‚2π‘π‘Ž2ξ€·+4π‘Ž1βˆ’π‘’2ξ€Έ2ξ‚€π‘Žπ‘Ÿξ‚(cos𝐸+cos𝑓)π‘π‘Žπ‘π‘’+ξ€·1βˆ’π‘’2ξ€Έ2(cos𝐸+cos𝑓)2𝑝𝑒2+ξ€·4π‘Ž1βˆ’π‘’2ξ€Έ2π‘’ξ‚€π‘Žπ‘Ÿξ‚ξ‚΅1sin𝑓1+×𝑝1+𝑒cosπ‘“π‘Žπ‘πœ”βˆ’ξ€·1βˆ’π‘’2ξ€Έ1/2π‘π‘Žπ‘π‘€ξ‚„+2ξ€·1βˆ’π‘’2ξ€Έ2𝑒1(cos𝐸+cos𝑓)1+𝑝1+𝑒cos𝑓×sinπ‘“π‘’π‘πœ”βˆ’βˆš1βˆ’π‘’2𝑝𝑒𝑝𝑀+1βˆ’π‘’2𝑒11+𝑝1+𝑒cos𝑓sinπ‘“πœ”βˆ’βˆš1βˆ’π‘’2𝑝𝑀2+12ξ‚€π‘Ÿπ‘Žξ‚2𝑝𝐼2+𝑝ΩsinπΌβˆ’π‘πœ”ξ‚cot𝐼2ξ‚Ή+12ξ‚€π‘Ÿπ‘Žξ‚2𝑝cos2(πœ”+𝑓)𝐼2βˆ’ξ‚€π‘Ξ©sinπΌβˆ’π‘πœ”ξ‚cot𝐼2ξ‚Ή+ξ‚€π‘Ÿπ‘Žξ‚2sin2(πœ”+𝑓)𝑝𝐼𝑝ΩsinπΌβˆ’π‘πœ”ξ‚ξ‚Ό,cot𝐼(3.9) with the disturbing function π‘‰ξ…ž2 given by π‘‰ξ…ž2=πœ‡π‘Žξ‚€π‘Žπ‘’π‘Žξ‚2𝐽2ξ‚»ξ‚€π‘Žπ‘Ÿξ‚312βˆ’34sin2𝐼+34ξ‚€π‘Žπ‘Ÿξ‚3sin2𝐼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 𝐻0 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, 𝐹0=π‘›ξ…žπ‘ξ…žπ‘€.(4.2)𝐹0 denotes the new undisturbed Hamiltonian. Now, consider the canonical system described by 𝐹0: π‘‘π‘Žξ…žπ‘‘π‘‘=0,π‘‘π‘’ξ…žπ‘‘π‘‘=0,π‘‘πΌξ…žπ‘‘π‘‘=0,π‘‘Ξ©ξ…žπ‘‘π‘‘=0,π‘‘πœ”ξ…žπ‘‘π‘‘=0,π‘‘π‘€ξ…žπ‘‘π‘‘=π‘›ξ…ž,π‘‘π‘ξ…žπ‘Ž=3𝑑𝑑2π‘›ξ…žπ‘Žξ…žπ‘ξ…žπ‘€,π‘‘π‘ξ…žπ‘’π‘‘π‘‘=0,π‘‘π‘ξ…žπΌπ‘‘π‘‘=0,π‘‘π‘ξ…žΞ©π‘‘π‘‘=0,π‘‘π‘ξ…žπœ”π‘‘π‘‘=0,π‘‘π‘ξ…žπ‘€π‘‘π‘‘=0,(4.3) general solution of which is given by π‘Žξ…ž=π‘Žξ…ž0,π‘’ξ…ž=π‘’ξ…ž0,πΌξ…ž=πΌξ…ž0,Ξ©ξ…ž=Ξ©ξ…ž0,πœ”ξ…ž=πœ”ξ…ž0,π‘€ξ…ž=π‘€ξ…ž0+π‘›ξ…žξ€·π‘‘βˆ’π‘‘0ξ€Έ,π‘ξ…žπ‘Ž=π‘ξ…žπ‘Ž0+32π‘›ξ…žξ€·π‘‘βˆ’π‘‘0ξ€Έπ‘Žξ…žπ‘ξ…žπ‘€,π‘ξ…žπ‘’=π‘ξ…žπ‘’0,π‘ξ…žπΌ=π‘ξ…žπΌ0,π‘ξ…žΞ©=π‘ξ…žΞ©0,π‘ξ…žπœ”=π‘ξ…žπœ”0,π‘ξ…žπ‘€=π‘ξ…žπ‘€0.(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 πœ•π‘†1πœ•π‘‘=𝐻𝐽2+π»π›Ύβˆ—βˆ’πΉ1.(4.5) According to [29], the mean value of 𝐻𝐽2+π»π›Ύβˆ— must be calculated, and, 𝑆1 is obtained through integration of the remaining part. 𝐹1 and 𝑆1 are given by the following equations: 𝐹13=βˆ’2π‘›ξ…žπ½2ξ‚€π‘Žπ‘’π‘Žξ…žξ‚2ξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’2cosπΌξ…žπ‘ξ…žΞ©+12ξ€·1βˆ’5cos2πΌξ…žξ€Έπ‘ξ…žπœ”ξ‚‡+π‘Žξ…žξƒ―2πœ‡4π‘Žξ…ž2π‘π‘Žξ…ž2+52ξ€·1βˆ’π‘’ξ…ž2ξ€Έπ‘π‘’ξ…ž2+ξ€·5βˆ’4π‘’ξ…ž2ξ€Έ2π‘’ξ…ž2π‘πœ”ξ…ž2+π‘πΌξ…ž22ξ€·1βˆ’π‘’ξ…ž2ξ€Έ31+2π‘’ξ…ž2+52π‘’ξ…ž2cos2πœ”ξ…žξ‚„+5π‘’ξ…ž2sin2πœ”ξ…ž2ξ€·1βˆ’π‘’ξ…ž2ξ€Έπ‘ξ…žπΌξ‚΅π‘ξ…žΞ©sinπΌξ…žβˆ’cotπΌξ…žπ‘ξ…žπœ”ξ‚Ά+12ξ€·1βˆ’π‘’ξ…ž2ξ€Έξ‚΅π‘ξ…žΞ©sinπΌξ…žβˆ’cotπΌξ…žπ‘ξ…žπœ”ξ‚Ά231+2π‘’ξ…ž2ξ‚βˆ’52π‘’ξ…ž2cos2πœ”ξ…žξ‚„ξƒ°,𝑆1=𝐽2ξ‚€π‘Žπ‘’π‘Žξ…žξ‚2ξ‚€3ξƒ―ξƒ―1βˆ’2sin2πΌξ…žξ‚ξƒ¬ξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Ά3βˆ’ξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’3/2ξƒ­+32sin2πΌξ…žξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Ά3ξ€·πœ”cos2ξ…ž+π‘“ξ…žξ€Έξƒ°π‘Žξ…žπ‘ξ…žπ‘Ž+ξƒ―32sin2πΌξ…žπ‘’ξ…žξ€·1βˆ’π‘’ξ…ž2ξ€Έξ‚Έβˆ’12ξ€·πœ”cos2ξ…ž+π‘“ξ…žξ€Έβˆ’π‘’ξ…ž2ξ‚€ξ€·cos2πœ”ξ…ž+π‘“ξ…žξ€Έ+13ξ€·cos2πœ”ξ…ž+3π‘“ξ…žξ€Έξ‚ξ‚Ή+ξ€·1βˆ’π‘’ξ…ž2ξ€Έπ‘’ξ…ž12βˆ’34sin2πΌξ…žξ‚ξƒ¬ξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Ά3βˆ’ξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’3/2ξƒ­+34sin2πΌξ…žξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Ά3ξ€·πœ”cos2ξ…ž+π‘“ξ…žξ€Έξƒ°π‘ξ…žπ‘’+ξƒ―34sin2πΌξ…žξ€·1βˆ’π‘’ξ…ž2ξ€Έ2ξ‚Έ12ξ€·πœ”cos2ξ…ž+π‘“ξ…žξ€Έ+π‘’ξ…ž2ξ‚€ξ€·cos2πœ”ξ…ž+π‘“ξ…žξ€Έ+13ξ€·cos2πœ”ξ…ž+3π‘“ξ…žξ€Έξ‚ξ‚Ήξƒ°π‘ξ…žπΌ+ξƒ―32cosπΌξ…žξ€·1βˆ’π‘’ξ…ž2ξ€Έ2ξ‚ƒβˆ’ξ€·π‘“ξ…žβˆ’π‘€ξ…ž+π‘’ξ…žsinπ‘“ξ…žξ€Έ+12ξ€·πœ”sin2ξ…ž+π‘“ξ…žξ€Έ+π‘’ξ…ž2ξ‚€ξ€·sin2πœ”ξ…ž+π‘“ξ…žξ€Έ+13ξ€·sin2πœ”ξ…ž+3π‘“ξ…žξ€Έξ‚ξ‚Ήξƒ°π‘ξ…žΞ©+ξƒ―34ξ€·5cos2πΌξ…žξ€Έβˆ’1ξ€·1βˆ’π‘’ξ…ž2ξ€Έ2ξ€·π‘“ξ…žβˆ’π‘€ξ…ž+π‘’ξ…žsinπ‘“ξ…žξ€Έ+14ξ€·3cos2πΌξ…žξ€Έβˆ’1π‘’ξ…žξ€·1βˆ’π‘’ξ…ž2ξ€Έξƒ¬ξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Ά2ξ€·1βˆ’π‘’ξ…ž2ξ€Έ+ξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Άξƒ­+1sinπ‘“ξ…ž+38sin2πΌξ…žπ‘’ξ…žξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’ξ‚΅π‘Žξƒ¬ξƒ¬ξ…žπ‘Ÿξ…žξ‚Ά2ξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’ξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Άξƒ­ξ€·+1sin2πœ”ξ…ž+π‘“ξ…žξ€Έ+ξƒ¬ξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Ά2ξ€·1βˆ’π‘’ξ…ž2ξ€Έ+ξ‚΅π‘Žξ…žπ‘Ÿξ…žξ‚Ά+13ξƒ­ξ€·sin2πœ”ξ…ž+3π‘“ξ…žξ€Έξƒ­+38ξ€·3βˆ’5cos2πΌξ…žξ€Έξ€·1βˆ’π‘’ξ…ž2ξ€Έ2ξ‚ƒξ€·πœ”sin2ξ…ž+π‘“ξ…žξ€Έ+π‘’ξ…žξ‚€ξ€·sin2πœ”ξ…ž+π‘“ξ…žξ€Έ+13ξ€·sin2πœ”ξ…ž+3π‘“ξ…žξ€Έπ‘ξ‚ξ‚„ξ‚‡ξ…žπœ”ξ‚‡+12ξƒŽπ‘Žξ…ž5πœ‡3ξƒ―8π‘’ξ…žsinπΈξ…žπ‘Žξ…ž2π‘π‘Žξ…ž2ξ€·+81βˆ’π‘’ξ…ž2ξ€ΈsinπΈξ…žπ‘Žξ…žπ‘ξ…žπ‘Žπ‘ξ…žπ‘’βˆ’8√1βˆ’π‘’ξ…ž2π‘’ξ…žcosπΈξ…žπ‘ξ…žπ‘Žπ‘ξ…žπœ”+ξ€·1βˆ’π‘’ξ…ž2ξ€Έξ‚ƒβˆ’54π‘’ξ…žsinπΈξ…ž+34sin2πΈξ…žβˆ’1𝑒12ξ…žsin3πΈξ…žξ‚„π‘π‘’ξ…ž2+√1βˆ’π‘’ξ…ž2π‘’ξ…žξ‚ƒ52π‘’ξ…žcosπΈξ…žβˆ’12ξ€·3βˆ’π‘’ξ…ž2ξ€Έcos2πΈξ…ž+16π‘’ξ…žcos3πΈξ…žξ‚„π‘ξ…žπ‘’π‘ξ…žπœ”+ξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’1ξƒ¬π‘πΌξ…ž2+ξ‚΅π‘ξ…žΞ©sinπΌξ…žβˆ’π‘ξ…žπœ”cotπΌξ…žξ‚Ά2ξƒ­Γ—ξ‚ƒξ‚€βˆ’π‘’ξ…ž+38π‘’ξ…ž3sinπΈξ…ž+38π‘’ξ…ž2sin2πΈξ…žβˆ’1𝑒24ξ…ž3sin3πΈξ…žξ‚„+ξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’1Γ—ξƒ¬π‘πΌξ…ž2cos2πœ”ξ…ž+2π‘ξ…žπΌξ‚΅π‘ξ…žΞ©sinπΌξ…žβˆ’π‘ξ…žπœ”cotπΌξ…žξ‚Άsin2πœ”ξ…žβˆ’ξ‚΅π‘ξ…žΞ©sinπΌξ…žβˆ’π‘ξ…žπœ”cotπΌξ…žξ‚Ά2cos2πœ”ξ…žξƒ­Γ—βˆ’54π‘’ξ…ž+58π‘’ξ…ž3sinπΈξ…ž+ξ‚€14+18π‘’ξ…ž2sin2πΈξ…ž+ξ‚€βˆ’1𝑒12ξ…ž+1𝑒24ξ…ž3sin3πΈξ…žξ‚„+ξ€·1βˆ’π‘’ξ…ž2ξ€Έβˆ’1/2ξ‚Έβˆ’π‘πΌξ…ž2sin2πœ”ξ…ž+2π‘ξ…žπΌξ‚΅π‘ξ…žΞ©sinπΌξ…žβˆ’π‘ξ…žπœ”cotπΌξ…žξ‚Άcos2πœ”ξ…ž+ξ‚΅π‘ξ…žΞ©sinπΌξ…žβˆ’π‘ξ…žπœ”cotπΌξ…žξ‚Ά2sin2πœ”ξ…žξƒ­Γ—ξ‚ƒ54π‘’ξ…žcosπΈξ…žβˆ’ξ‚€14+14π‘’ξ…ž2cos2πΈξ…ž+1𝑒12ξ…žcos3πΈξ…žξ‚„+π‘πœ”ξ…ž2π‘’ξ…ž254π‘’ξ…žβˆ’π‘’ξ…ž3sinπΈξ…ž+ξ‚€βˆ’34+12π‘’ξ…ž2sin2πΈξ…ž+1𝑒12ξ…ž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 𝐹13=βˆ’2𝑛𝐽2ξ‚€π‘Žπ‘’π‘Žξ‚2cos𝐼𝑝Ω+π‘Žξƒ―2πœ‡4π‘Ž2𝑝2π‘Ž+52𝑝2β„Ž+𝑝2π‘˜ξ€Έ+12𝑝2𝐼+𝑝2Ξ©sin2𝐼,𝑆ξƒͺξƒ°1=32𝐽2ξ‚€π‘Žπ‘’π‘Žξ‚2ξ‚†π‘Žπ‘π‘Žsin2𝐼cos2β„“+π‘β„Žξ‚ƒ1cosβ„“+4sin2𝐼7βˆ’5cosβ„“+3cos3ℓ+π‘π‘˜ξ‚ƒ1sinβ„“+4sin2𝐼7βˆ’7sinβ„“+3+1sin3ℓ4𝑝𝐼1sin2𝐼cos2β„“+2𝑝Ω+1cos𝐼sin2β„“2ξƒŽπ‘Ž5πœ‡38π‘Žπ‘π‘Žξ€Ίπ‘β„Žsinβ„“βˆ’π‘π‘˜ξ€»βˆ’1cosβ„“23π‘β„Žπ‘π‘˜+𝑝Ω𝑝sin𝐼𝐼+1cos2β„“43𝑝2β„Žβˆ’π‘2π‘˜ξ€Έ+𝑝2πΌβˆ’π‘2Ξ©sin2𝐼.sin2β„“(5.3) Double prime is omitted to simplify the notation.

For transfers between close orbits, 𝐹1 and 𝑆1 can be linearized around a suitable reference orbit, and an approximate first-order analytical solution can be determined. This solution is given by Ξ”π‘₯=𝐴𝑝0+𝐡,(5.4) where Ξ”π‘₯=[Ξ”π›ΌΞ”β„ŽΞ”π‘˜Ξ”πΌΞ”Ξ©ξ…ž]𝑇,𝑝0 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π‘˜π›Όπ‘Žπ‘˜β„Žπ‘Žπ‘˜π‘˜00000π‘ŽπΌπΌπ‘ŽπΌΞ©β€²000π‘ŽΞ©β€²πΌπ‘ŽΞ©β€²Ξ©β€²βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ¦βŽ‘βŽ’βŽ’βŽ’βŽ£π‘,(5.5)𝐡=π›Όπ‘β„Žπ‘π‘˜π‘πΌπ‘Ξ©β€²βŽ€βŽ₯βŽ₯βŽ₯⎦,(5.6) where π‘Žπ›Όπ›ΌξƒŽ=4π‘Ž5πœ‡3β„“|||||ℓ𝑓ℓ0,π‘Žπ›Όβ„Ž=π‘Žβ„Žπ›ΌξƒŽ=4π‘Ž5πœ‡3sinβ„“|||||ℓ𝑓ℓ0,π‘Žπ›Όπ‘˜=π‘Žπ‘˜π›ΌξƒŽ=βˆ’4π‘Ž5πœ‡3cosβ„“|||||ℓ𝑓ℓ0,π‘Žβ„Žβ„Ž=ξƒŽπ‘Ž5πœ‡3523β„“+4sin2β„“ξ‚„|||||ℓ𝑓ℓ0,π‘Žβ„Žπ‘˜=π‘Žπ‘˜β„Ž3=βˆ’4ξƒŽπ‘Ž5πœ‡3cos2β„“|||||ℓ𝑓ℓ0,π‘Žπ‘˜π‘˜=ξƒŽπ‘Ž5πœ‡3523β„“βˆ’4sin2β„“ξ‚„|||||ℓ𝑓ℓ0,π‘ŽπΌπΌ=ξƒŽπ‘Ž5πœ‡3121β„“+4sin2β„“ξ‚„|||||ℓ𝑓ℓ0,π‘ŽπΌΞ©β€²=π‘ŽΞ©β€²πΌ1=βˆ’4ξƒŽπ‘Ž5πœ‡3cos2β„“|||||ℓ𝑓ℓ0,π‘ŽΞ©β€²Ξ©β€²=ξƒŽπ‘Ž5πœ‡3121β„“βˆ’4sin2β„“ξ‚„|||||ℓ𝑓ℓ0,𝑏𝛼=πœ€sin2𝐼cos2β„“|||ℓ𝑓ℓ0,π‘β„Žξ‚ƒ=πœ€cos1β„“+4sin2πΌξ‚€βˆ’5cos7β„“+3cos3β„“|||ℓ𝑓ℓ0,π‘π‘˜ξ‚ƒ=πœ€sin1β„“+4sin2πΌξ‚€βˆ’7sin7β„“+3sin3β„“|||ℓ𝑓ℓ0,𝑏𝐼=πœ€4sin2𝐼cos2β„“|||ℓ𝑓ℓ0,π‘Ξ©β€²πœ€=βˆ’2sin2𝐼1β„“βˆ’2sin2ℓ|||ℓ𝑓ℓ0,(5.7) with ℓ𝑓=β„“0+𝑛(π‘‘π‘“βˆ’π‘‘0),𝑑0 is the initial time, 𝑑𝑓 is the final time, and πœ€=(3/2)𝐽2(π‘Žπ‘’/π‘Ž)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 𝐽2 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=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 1Δ𝐽=2ξ€½π‘Žπ›Όπ›Όπ‘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ξƒŽπœ‡3π‘Ž5Δ𝛼Δℓ,π‘β„Ž=25ξƒŽπœ‡3π‘Ž5ξ‚΅Ξ”β„ŽΞ”β„“ξ‚Ά,π‘π‘˜=25ξƒŽπœ‡3π‘Ž5ξ‚΅Ξ”π‘˜Ξ”β„“ξ‚Ά,π‘πΌξƒŽ=2πœ‡3π‘Ž5Δ𝐼Δℓ,π‘Ξ©β€²ξƒŽ=2πœ‡3π‘Ž5ΔΩsin𝐼Δℓ+πœ€2sin2𝐼.(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: πœΈβˆ—=πœ‡π‘Ž2ξ‚»25ξ‚΅Ξ”β„ŽΞ”β„“sinβ„“βˆ’Ξ”π‘˜Ξ”β„“cosβ„“ξ‚Άπžπ‘Ÿ+ξ‚΅12Δ𝛼Δℓ+45ξ‚΅Ξ”β„ŽΞ”β„“cosβ„“+Ξ”π‘˜Ξ”β„“sinβ„“πžξ‚Άξ‚Άπ‘ ξ‚΅+2Δ𝐼Δℓcosξ‚Έβ„“+ΔΩΔℓ+πœ€cos𝐼sin𝐼sinβ„“ξ‚Άπžπ‘€ξ‚Ό,(5.11)Δ𝐽Δℓ=14ξƒŽπœ‡3π‘Ž5ξƒ―12Δ𝛼Δℓ2+45ξƒ©Ξ”β„Ž2+Ξ”π‘˜2Ξ”β„“2ξƒͺ+4Δ𝐼Δℓ2+ΔΩΔℓ+πœ€cos𝐼2sin2𝐼.ξƒͺξƒ°(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 𝐹1=32𝑛𝐽2ξ‚€π‘Žπ‘’π‘Žξ‚2ξ€·π‘„π‘π‘ƒβˆ’π‘ƒπ‘π‘„ξ€Έ+π‘Žξ‚†2πœ‡4π‘Ž2𝑝2π‘Ž+52𝑝2πœ‰+𝑝2πœ‚ξ‚+18𝑝2𝑃+𝑝2𝑄,𝑆(5.15)1=32𝐽2ξ‚€π‘Žπ‘’π‘Žξ‚2ξ‚†π‘πœ‰cosπœ†+π‘πœ‚1sinπœ†+2𝑄𝑝𝑃+𝑃𝑝𝑄1sin2πœ†+2ξ€·π‘ƒπ‘π‘ƒβˆ’π‘„π‘π‘„ξ€Έξ‚‡+1cos2πœ†2ξƒŽπ‘Ž5πœ‡38π‘Žπ‘π‘Žξ€·π‘πœ‰sinπœ†βˆ’π‘πœ‚ξ€Έ+3cosπœ†4𝑝2πœ‰βˆ’π‘2πœ‚ξ‚+1𝑝162π‘ƒβˆ’π‘2π‘„βˆ’ξ‚ƒ3sin2πœ†2π‘πœ‰π‘πœ‚+18𝑝𝑃𝑝𝑄.cos2πœ†(5.16) Double prime is omitted to simplify the notation.

For transfers between close orbits, 𝐹1 and 𝑆1 can be linearized around a suitable reference orbit, and an approximate first-order analytical solution can be determined. This solution is given by Ξ”π‘₯=𝐢𝑝0+𝐷,(5.17) where Ξ”π‘₯=[Ξ”π›ΌΞ”πœ‰Ξ”πœ‚Ξ”π‘ƒΞ”π‘„]𝑇,𝑝0 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πœ‚π›Όπ‘πœ‚πœ‰π‘πœ‚πœ‚00000𝑐𝑃𝑃𝑐𝑃𝑄000π‘π‘„π‘ƒπ‘π‘„π‘„βŽ€βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯⎦⎑⎒⎒⎒⎣0𝑑,(5.18)𝐷=πœ‰π‘‘πœ‚π‘‘π‘ƒπ‘‘π‘„βŽ€βŽ₯βŽ₯βŽ₯⎦,(5.19) where π‘π›Όπ›ΌξƒŽ=4π‘Ž5πœ‡3πœ†|||||πœ†π‘“πœ†0,π‘π›Όπœ‰=π‘πœ‰π›ΌξƒŽ=4π‘Ž5πœ‡3sinπœ†|||||πœ†π‘“πœ†0,π‘π›Όπœ‚=π‘πœ‚π›ΌξƒŽ=βˆ’4π‘Ž5πœ‡3cosπœ†|||||πœ†π‘“πœ†0,π‘πœ‰πœ‰=ξƒŽπ‘Ž5πœ‡3523πœ†+4sin2πœ†ξ‚„|||||πœ†π‘“πœ†0,π‘πœ‰πœ‚=π‘πœ‚πœ‰3=βˆ’4ξƒŽπ‘Ž5πœ‡3cos2πœ†|||||πœ†π‘“πœ†0,π‘πœ‚πœ‚=ξƒŽπ‘Ž5πœ‡3523πœ†βˆ’4sin2πœ†ξ‚„|||||πœ†π‘“πœ†0,𝑐𝑃𝑃=18ξƒŽπ‘Ž5πœ‡3ξ‚ƒΞ”πœ†cosπœ€Ξ”1πœ†+2ξ‚€sin(2+πœ€)πœ†π‘“βˆ’πœ€πœ†0ξ‚βˆ’12ξ‚€sin(2+πœ€)πœ†0βˆ’πœ€πœ†π‘“,𝑐𝑃𝑄=18ξƒŽπ‘Ž5πœ‡3ξ‚ƒΞ”πœ†sinπœ€Ξ”1πœ†βˆ’2ξ‚€cos(2+πœ€)πœ†π‘“βˆ’πœ€πœ†0+12ξ‚€cos(2+πœ€)πœ†0βˆ’πœ€πœ†π‘“,𝑐𝑄𝑃=18ξƒŽπ‘Ž5πœ‡3ξ‚ƒβˆ’Ξ”πœ†sinπœ€Ξ”1πœ†βˆ’2ξ‚€cos(2+πœ€)πœ†π‘“βˆ’πœ€πœ†0+12ξ‚€cos(2+πœ€)πœ†0βˆ’πœ€πœ†π‘“,𝑐𝑄𝑄=18ξƒŽπ‘Ž5πœ‡3ξ‚ƒΞ”πœ†cosπœ€Ξ”1πœ†βˆ’2ξ‚€(sin2+πœ€)πœ†π‘“βˆ’πœ€πœ†0+12ξ‚€(sin2+πœ€)πœ†0βˆ’πœ€πœ†π‘“,π‘‘ξ‚ξ‚„πœ‰=32𝐽2ξ‚€π‘Žπ‘’π‘Žξ‚2cosπœ†||||πœ†π‘“πœ†0,π‘‘πœ‚=32𝐽2ξ‚€π‘Žπ‘’π‘Žξ‚2sinπœ†||||πœ†π‘“πœ†0,𝑑𝑃=𝑃0cosπœ€Ξ”1πœ†βˆ’2ξ‚€πœ€πœ€cosπœ†π‘“βˆ’(2+πœ€)πœ†0+12ξ‚€πœ€πœ€cosπœ†0βˆ’(2+πœ€)πœ†π‘“ξ‚‡ξ‚ξ‚„βˆ’1+𝑄0sinπœ€Ξ”1πœ†+2ξ‚€πœ€πœ€sinπœ†π‘“βˆ’(2+πœ€)πœ†0ξ‚βˆ’12ξ‚€πœ€πœ€sinπœ†0βˆ’(2+πœ€)πœ†π‘“,𝑑𝑄=𝑄0cosπœ€Ξ”1πœ†+2ξ‚€πœ€πœ€cosπœ†π‘“βˆ’(2+πœ€)πœ†0ξ‚βˆ’12ξ‚€πœ€πœ€cosπœ†0βˆ’(2+πœ€)πœ†π‘“ξ‚‡ξ‚ξ‚„βˆ’1+𝑃0ξ‚ƒβˆ’sinπœ€Ξ”1πœ†+2ξ‚€πœ€πœ€sinπœ†π‘“βˆ’(2+πœ€)πœ†0ξ‚βˆ’12ξ‚€πœ€πœ€sinπœ†0βˆ’(2+πœ€)πœ†π‘“.(5.20) The adjoint variables are π‘π‘Ž=π‘π‘Ž0,π‘πœ‰=π‘πœ‰0,π‘πœ‚=π‘πœ‚0,𝑝𝑃=𝑝𝑃0cosπœ€Ξ”πœ†+𝑝𝑄0sinπœ€Ξ”π‘πœ†,𝑄=βˆ’π‘π‘ƒ0sinπœ€Ξ”πœ†+𝑝𝑄0cosπœ€Ξ”πœ†,(5.21) where Ξ”πœ†=πœ†π‘“βˆ’πœ†0=𝑛(π‘‘π‘“βˆ’π‘‘0). 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 𝐽2 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=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𝑝𝑃cosπœ†+π‘π‘„ξ€Έπžsinπœ†π‘€ξ‚‡,1Δ𝐽=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ξƒŽπœ‡3π‘Ž5ξ‚΅Ξ”π›ΌΞ”πœ†ξ‚Ά,π‘πœ‰=25ξƒŽπœ‡3π‘Ž5ξ‚΅Ξ”πœ‰Ξ”πœ†ξ‚Ά,π‘πœ‚=25ξƒŽπœ‡3π‘Ž5ξ‚΅Ξ”πœ‚Ξ”πœ†ξ‚Ά,𝑝𝑃=ξƒŽπœ‡3π‘Ž58Ξ”πœ†ξ‚ƒξ€·π‘ƒ0ξ€Έξ‚€+Δ𝑃cosπœ€Ξ”πœ†ξ‚βˆ’ξ€·π‘„0ξ€Έξ‚€+Δ𝑄sinπœ€Ξ”πœ†ξ‚βˆ’π‘ƒ0ξ‚„,𝑝𝑄=ξƒŽπœ‡3π‘Ž58Ξ”πœ†ξ‚ƒξ€·π‘„0ξ€Έξ‚€+Δ𝑄cosπœ€Ξ”πœ†ξ‚+𝑃0ξ€Έξ‚€+Δ𝑃sinπœ€Ξ”πœ†ξ‚βˆ’π‘„0ξ‚„.(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: πœΈβˆ—=πœ‡π‘Ž2ξ‚»25ξ‚΅Ξ”πœ‰Ξ”πœ†sinπœ†βˆ’Ξ”πœ‚Ξ”πœ†cosπœ†ξ‚Άπžπ‘Ÿ+ξ‚΅12Ξ”π›ΌΞ”πœ†+45ξ‚΅Ξ”πœ‰Ξ”πœ†cosπœ†+Ξ”πœ‚Ξ”πœ†sinπœ†πžξ‚Άξ‚Άπ‘ +4Ξ”πœ†ξ‚ƒπ‘ƒπ‘“ξ‚€cosπœ†+πœƒπ‘“ξ‚βˆ’π‘ƒ0ξ‚€cosπœ†+πœƒ0+𝑄𝑓sinπœ†+πœƒπ‘“ξ‚βˆ’π‘„0ξ‚€sinπœ†+πœƒ0πžξ‚ξ‚„π‘€ξ‚Ό,Ξ”π½Ξ”πœ†=14ξƒŽπœ‡3π‘Ž5ξƒ―12ξ‚΅Ξ”π›ΌΞ”πœ†ξ‚Ά2+45ξƒ©Ξ”πœ‰2+Ξ”πœ‚2Ξ”πœ†2ξƒͺξ‚€+16cosπœ€Ξ”πœ†ξ‚ξƒ―ξƒ©Ξ”π‘ƒ2+Δ𝑄2Ξ”πœ†2ξƒͺ+4Ξ”πœ†2ξ‚΅sinπœ€Ξ”πœ†2×𝑃𝑓𝑃0ξ‚΅sinπœ€Ξ”πœ†2ξ‚Άβˆ’π‘„0ξ‚΅cosπœ€Ξ”πœ†2ξ‚Άξ‚Ά+𝑄𝑓𝑃0ξ‚΅cosπœ€Ξ”πœ†2ξ‚Ά+𝑄0ξ‚΅sinπœ€Ξ”πœ†2,ξ‚Άξ‚Άξ‚Ήξ‚Όξ‚Ό(5.24) with πœƒ0=πœ€π‘›(π‘‘βˆ’π‘‘0) 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), 𝐹1 is written as 𝐹13=βˆ’2𝑛𝐽2ξ‚€π‘Žπ‘’π‘Žξ‚2ξ€·1βˆ’πœ‰2βˆ’πœ‚2ξ€Έβˆ’2𝑃1βˆ’22+𝑄2ξ€Έξ€Έξ€Ίξ€·πœ‰π‘πœ‚βˆ’πœ‚π‘πœ‰ξ€Έ+ξ€·π‘ƒπ‘π‘„βˆ’π‘„π‘π‘ƒ+1ξ€Έξ€»2𝑃1βˆ’51βˆ’22+𝑄2ξ€Έξ€Έ2ξ‚ξ€·πœ‰π‘πœ‚βˆ’πœ‚π‘πœ‰ξ€Έξ‚‡+π‘Žξ‚†2πœ‡4π‘Ž2𝑝2π‘Ž+𝑝2πœ‰ξ‚€52βˆ’52πœ‰2βˆ’2πœ‚2+𝑝2πœ‚ξ‚€52βˆ’52πœ‚2βˆ’2πœ‰