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

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.

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/𝑔₀, is between 10⁻⁴ and 10⁻². For such system, the fuel consumption is described by the variable 𝐽 defined as 𝐽=12^𝑡𝑓_𝑡0𝛾²𝑑𝑡,(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 the 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 𝑚_𝑓.

Let us consider the motion of a space vehicle 𝑀, powered by a limited-power engine in a Newtonian central gravity field. At time 𝑡, the state of the vehicle is defined by the position vector 𝐫(𝑡), the velocity vector 𝐯(𝑡) and the consumption variable 𝐽(𝑡). 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 𝑀 from the initial state (𝐫₀,𝐯₀,0) at time 𝑡₀ to the final state (𝐫_𝑓,𝐯_𝑓,𝐽_𝑓) at time 𝑡_𝑓, such that the final consumption variable 𝐽_𝑓 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 that 𝑡_𝑓−𝑡₀ and the position of the vehicle in the initial orbit are specified. The state equations are 𝑑𝐫𝑑𝑑𝐯𝑡=𝐯,𝑑𝑡=−𝑟𝜇³𝑑𝐽𝐫+𝜸,𝑑𝑡=12𝛾²,(2.3) 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 𝐻 reaches its maximum. The Hamiltonian 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, it follows that 𝜸^∗=−𝐩^𝑣𝑝_𝐽.(2.5) The optimal thrust acceleration 𝜸^∗ is modulated [1]. Then, the optimal trajectories are governed by the maximum Hamiltonian 𝐻^∗, 𝐻^∗=𝐩_𝑟•𝐯−𝐩_𝑣•𝜇𝑟³𝐫−𝑝^2𝑣2𝑝_𝐽.(2.6)Since 𝐽 is an ignorable variable, it is determined by simple integration, and its adjoint 𝑝_𝐽 is a first integral, whose value is obtained from the transversality conditions: 𝑝_𝐽=−1. Thus, 𝐻^∗=𝐩_𝑟•𝐯−𝐩_𝑣•𝜇𝑟³𝐫+𝑝^2𝑣2.(2.7)

This Hamiltonian can be decomposed into two parts: 𝐻₀=𝐩_𝑟•𝐯−𝐩_𝑣•𝜇𝑟³𝐻𝐫,_𝛾∗=𝑝^2𝑣2,(2.8)where 𝐻₀ is the undisturbed Hamiltonian and 𝐻_𝛾∗ is the disturbing function concerning the optimal thrust acceleration. The undisturbed Hamiltonian 𝐻₀ 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 Hamiltonian 𝐻₀, that is, the problem defined by the differential equations for the state variables 𝐫 and 𝐯, which describe the classical two-body problem, and the differential equations for the adjoint variables 𝐩_𝑟 and 𝐩_𝑣. These equations are given by 𝑑𝐫𝑑𝑡=𝐯𝑑𝑑𝐩𝑡^𝐫=𝑟𝜇³𝐩^𝑣𝑑𝐯−3𝐩_𝑣•𝐞_𝐫𝐞^𝐫,𝑑𝑡=−𝑟𝜇³𝐫𝑑𝐩^𝑣𝑑𝑡=−𝐩^𝑟,(3.1) where 𝐞_𝐫 is the unit vector pointing radially outward of the moving frame of reference.

The general solution of the differential equations for the state variables 𝐫 and 𝐯 is well-known from the classical two-body problem [8]. In the coplanar case, it is given by 𝐫=𝑎1−𝑒²1+𝑒cos𝑓𝐞^𝐫𝐯=⎷,(3.2)𝜇𝑎1−𝑒²𝑒sin𝑓𝐞^𝐫+1+𝑒cos𝑓𝐞_𝐬,(3.3)𝐞_𝐬where 𝑎 is the unit vector along circumferential direction of the moving frame of reference; 𝑒 is the semimajor axis, 𝑓 is the eccentricity, and 𝐞_𝐫 is the true anomaly. The unit vectors 𝐞_𝐬 and 𝐞𝐞_𝐫=cos(𝜔+𝑓)𝐢+sin(𝜔+𝑓)𝐣,_𝐬=−sin(𝜔+𝑓)𝐢+cos(𝜔+𝑓)𝐣(3.4) are written as function of the orbital elements as follows: 𝐢 where 𝐣 and 𝜔 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 variables 𝐩_𝑣 and 𝑎=𝑟𝑒2−𝑟𝑣²/𝜇,(3.5)₂=1−ℎ²cos𝜔=𝐢𝜇•𝑎𝐞,(3.6)cos𝐸=1𝑒,(3.7)𝑒1−𝑎𝑟,(3.8) 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): 𝐞𝐡𝐫𝐯where the eccentricity vector 𝐞=1𝜇𝑣²−𝜇𝐡=𝑟𝐫×𝐫𝐯−.𝐫•𝐯𝐯,(3.9) and the angular momentum vector ×, shown in Figure 2, are given, as function of the Cartesian elements 𝑓 and 𝐸, by the following equations: tan𝑓2=1+𝑒1−𝑒tan𝐸2.(3.10) Here, the symbol 𝐞 denotes the cross product. The true anomaly 𝐡 has been replaced by eccentric anomaly 𝛿𝐫. These anomalies are related through the equation 𝛿𝐯

Figure 2 

The vectors 𝛿𝛿𝐫𝐯==𝛿𝛿𝑢𝜉𝐞𝐞𝐫𝐫+𝛿𝜂𝐞𝐬,+𝛿𝑣𝐞𝐬.(3.11) and 𝑎.

Now, consider the variations in the Cartesian elements, 𝑒 and 𝜔, given in the moving frame of reference by 𝐸

The variations of the orbital elements—𝛿𝐫, 𝛿𝐯, 𝛿𝑎=2𝑎𝑟²𝐫•𝛿𝐫𝛿𝑒=−𝐡•𝛿𝐡𝑟+2𝜇𝑎²𝐯•𝛿𝐯,𝛿𝜔=𝛿𝑒𝜇𝑎𝑒+2𝜇ℎ𝑎²²𝑒𝛿𝑎,𝛿𝐸=𝑒1cot𝜔−𝑒𝐢s•in𝛿𝐞𝜔,𝑒sin𝐸𝐫•𝑟𝑎𝛿𝐫−𝑎𝑟²𝛿𝑎+cos𝐸𝛿𝑒,(3.12), and 𝐞—induced by the variations in the Cartesian elements, 𝐡 and 𝜇𝛿𝐞=𝑣^2𝑠𝛿𝜉−𝑣_𝑠𝑣_𝑟𝛿𝜂+2𝑟𝑣_𝑠𝛿𝑣𝐞_𝐫+−𝑣_𝑠𝑣_𝑟𝛿𝜉+𝑣^2𝑟−𝜇𝑟𝛿𝜂−𝑟𝑣^𝑠𝛿𝐡=𝑣_𝑠𝛿𝑢−𝑟𝑣_𝑟𝛿𝑣𝐞_𝐬,(3.13)𝛿𝜉−𝑣_𝑟𝛿𝜂+𝑟𝛿𝑣𝐞_𝐰,(3.14), are obtained straightforwardly from (3.5) through (3.8) and are given by 𝐞_𝐰 where the variations of the vectors 𝑣_𝑟 and 𝑣_𝑠 are written as 𝑎𝑒being the unit vector 𝜔 normal to the plane of the orbit. Here, 𝐸 and 𝛿𝐫 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—𝛿𝐯, 𝐸, 𝑀, and 𝑀=𝐸−𝑒sin𝐸.(3.15)—induced by the variations in the Cartesian elements, 𝛿𝑀=𝑟𝑎𝛿𝐸−sin𝐸𝛿𝑒.(3.16) and 𝛿𝑎=2𝑎𝑟²𝛿𝜉+2𝑒sin𝑓𝑛√1−𝑒²𝛿𝑢+2√1−𝑒²𝛿𝑒=𝑎1−𝑒₂𝑛𝑎𝑟𝛿𝑣,𝑟²cos𝐸𝛿𝜉+sin𝑓𝑎𝛿𝜂+√1−𝑒²𝑛𝑎sin𝑓𝛿𝑢+√1−𝑒²𝛿𝜔=sin𝑓𝑛𝑎cos𝐸+cos𝑓𝛿𝑣,𝑒𝑟𝛿𝜉−𝑎𝑒𝑒+1c−os𝑒𝑓₂𝛿𝜂−√1−𝑒²𝑛𝑎𝑒cos𝑓𝛿𝑢+√1−𝑒²𝛿𝑀=−1−𝑒₃𝑛𝑎𝑒sin𝑓1+1+𝑒1cos𝑓𝛿𝑣,cos𝐸𝑒𝑟√1−𝑒²sin𝑓𝛿𝜉+√1−𝑒²+1−𝑒₂𝑎𝑒cos𝑓𝛿𝜂𝑛𝑎𝑒cos𝑓−1+𝑒2c𝑒os𝑓𝛿𝑢−1−𝑒²𝑛𝑎𝑒sin𝑓1+1+𝑒1cos𝑓𝛿𝑣,(3.17). At this point, it is useful to replace the variation in the eccentric anomaly 𝑛=√𝜇/𝑎³ by the variation in the mean anomaly ⎡⎢⎢⎢⎢⎢⎣𝛿𝛿𝛿𝛿𝑀𝜔𝑎𝑒⎤⎥⎥⎥⎥⎥⎦=Δ₋₁⎡⎢⎢⎢⎢⎢⎣𝛿𝛿𝛿𝛿𝑣𝜂𝜉𝑢⎤⎥⎥⎥⎥⎥⎦,(3.18), obtained from the well-known Kepler equation Δ⁻¹This variation is given by 𝐩_𝐫Thus, 𝐩_𝐯 where 𝐩𝐩^𝐫_𝐯=Δ⁻¹^𝑇⎡⎢⎢⎢⎢⎢⎣𝑝^𝑎𝑝_𝑒𝑝_𝜔𝑝_𝑀⎤⎥⎥⎥⎥⎥⎦,(3.19) is the mean motion.

Equation (3.17) can be put in the form 𝐩_𝐫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 [7], the general solution of the differential equations for the adjoint variables 𝐩_𝐫=𝑎𝑟²2𝑎𝑝^𝑎+1−𝑒²cos𝐸𝑝^𝑒+𝑟𝑎sin𝑒𝑓𝑝^𝜔−1−𝑒³cos𝐸√1−𝑒²𝑝_𝑀𝐞^𝐫+sin𝑓𝑎𝑝^𝑒−𝑒+cos𝑓𝑎𝑒1−𝑒²𝑝^𝜔+√1−𝑒²cos𝑓𝑎𝑒𝑝^𝑀𝐞^𝐬𝐩_𝐯,=1𝑛𝑎√1−𝑒²2𝑎𝑒sin𝑓𝑝^𝑎+1−𝑒²sin𝑓𝑝^𝑒−1−𝑒²cos𝑓𝑒𝑝^𝜔+1−𝑒²^3/2𝑒cos𝑓−1+𝑒2c𝑒os𝑓𝑝^𝑀𝐞^𝐫+2𝑎1−𝑒²𝑎𝑟𝑝^𝑎+1−𝑒²cos𝑓+cos𝐸𝑝^𝑒+1−𝑒²sin𝑓𝑒1+1+𝑒1cos𝑓(𝑝^𝜔−√1−𝑒²𝑝_𝑀)𝐞_𝐬.(3.20) and (𝑟/𝑎) is given by (𝑟/𝑎)sin𝑓with (𝐫,𝐯,𝐩_𝐫,𝐩_𝐯) and (𝑎,𝑒,𝜔,𝑀,𝑝_𝑎,𝑝_𝑒,𝑝_𝜔,𝑝_𝑀) expressed in the moving frame of reference. Thus, it follows from (3.17) through (3.19) that 𝐻^∗𝐻𝐻₀=𝑛𝑝_𝑀,(3.21)_𝛾∗=1×2𝑛²1𝑎²1−𝑒²2(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−𝑒²𝑝_𝑎𝑝_𝑀]²𝑒(cos𝐸+cos𝑓)1+1+𝑒1cos𝑓sin𝑓[𝑝^𝑒+1−𝑒₂𝑝_𝜔−√1−𝑒²𝑝_𝑒𝑝_𝑀]𝑒1+1+𝑒1cos𝑓sin𝑓[𝑝^𝜔−√1−𝑒²𝑝_𝑀]².(3.22), 𝐻^∗, 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 𝐻^∗ and the orbital ones 𝐻^∗. The Hamiltonian 𝐻₀ is invariant with respect to this canonical transformation. Thus, 𝐻_𝛾∗

The new Hamiltonian function 𝐻₀, 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 𝐻_𝛾∗ 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 function 𝑎,𝑒,𝜔,𝑀,𝑝^𝑎,𝑝_𝑒,𝑝_𝜔,𝑝_𝑀→𝑎^,𝑒^,𝜔^,𝑀^,𝑝^𝑎,𝑝^𝑒,𝑝^𝜔,𝑝^𝑀.(4.1) 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 functions 𝐹₀=𝑛^𝑝^𝑀.(4.2) and 𝐹₀ have different order of magnitude: 𝑑𝑎^𝑑𝑡=0,𝑑𝑑𝑝𝑡^𝑎=32𝑎𝑛^𝑑𝑒_𝑝^𝑀,𝑑𝑑𝜔𝑡_=0,𝑑𝑑𝑝𝑡^𝑒=0,𝑑𝑀𝑑𝑡_=0,𝑑𝑑𝑝𝑡^𝜔=0,𝑑𝑡=𝑛,𝑑𝑑𝑝^𝑀𝑡=0,(4.3) has zero order and 𝑝𝑎^_𝑎=𝑎^0,𝑒^=𝑒^0,𝜔^=𝜔^0,𝑀^=𝑀^0+𝑛^𝑡−𝑡⁰,=𝑝^𝑎₀+32𝑛^𝑡−𝑡⁰𝑎^𝑝^𝑀,𝑝^𝑒=𝑝^𝑒₀,𝑝^𝜔=𝑝^𝜔₀,𝑝^𝑀=𝑝^𝑀₀.(4.4) 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.5) Following the algorithm described in the appendix, at order 0 one gets 𝐻_𝛾∗

Now, consider the undisturbed canonical system described by 𝐹₁ : 𝐻_𝛾∗ general solution of which is given by 𝐹₁=𝑎^2𝜇4𝑎^2𝑝^2𝑎+521−𝑒^2𝑝^2𝑒+5−4𝑒^22𝑒^2𝑆𝑝^2𝜔,(4.6)₁=12𝑎^5𝜇³8𝑒^sin𝐸^𝑎^2𝑝^𝑎2+81−𝑒^2sin𝐸^𝑎^𝑝^𝑎𝑝^𝑒−81−𝑒^2^1/2𝑒^cos𝐸^𝑎^𝑝^𝑎𝑝^𝜔+1−𝑒^2−54𝑒^sin𝐸^+34sin2𝐸^−112𝑒^sin3𝐸^𝑝^𝑒2+21−𝑒^2^1/2𝑒^54𝑒^cos𝐸^+14𝑒^2−3cos2𝐸^+112𝑒^cos3𝐸^𝑝^𝑒𝑝^𝜔+1𝑒^254−𝑒^2𝑒^sin𝐸^−1232−𝑒^2sin2𝐸^+112𝑒^sin3𝐸^𝑝^𝜔2.(4.7)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 𝑝^𝑀The functions 𝐹^=𝐹₁ and 𝑆=𝑆₁ are written in terms of the constants of integration. According to (A.11), the mean value of 𝐹^ must be calculated. After lengthy calculations, one gets 𝑎^,𝑒^,𝜔^,𝑝^𝑎,𝑝^𝑒,𝑝^𝜔→𝑎^,𝜙,𝜔^,𝑝^𝑎,𝑝_𝜙,𝑝^𝜔,(4.8) Terms factored by 𝑎𝑒^_=𝑎^,𝑝^𝑎=𝑝^𝑎,=sin𝜙,𝑝^𝑒=𝑝^𝜙𝜔cos𝜙_,=𝜔^,𝑝^𝜔=𝑝^𝜔.(4.9) 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, 𝐹^ and 𝐹^=𝑎^2𝜇4𝑎^2𝑝^𝑎2+52𝑝^𝜙2+52csc²𝜙−2𝑝^𝜔2.(4.10).

The new Hamiltonian 𝐹^ 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: 𝑡defined by the following equations: 𝜔The Hamiltonian function 𝑝𝑝^𝜔=𝐶₁,(4.11)_2𝜙𝑎+𝑝^𝜔2csc²𝜙=𝐶²²,(4.12)_2𝜇4𝑎^2𝑝^𝑎2+52𝑝^𝜙2+52csc²𝜙−2𝑝^𝜔2=E.(4.13) is invariant with respect to this transformation. Thus, E

The canonical system governed by the new Hamiltonian function 𝐸 has three first integrals (𝐹^ and 𝑊 are ignorable variables), 𝐶₁,𝐶₂,E(Constant 𝑎^,𝜙,𝜔^,𝑝^𝑎,𝑝_𝜙,𝑝^𝜔→^𝑊𝐶₁,𝐶₂,E,𝑝_𝐶1,𝑝_𝐶2,𝑝_E.(4.14) should not be confused with the eccentric anomaly 𝐹^.) These first integrals play important role in determining the general solution of the canonical system governed by the new Hamiltonian function 𝑝^𝑎=𝜕𝑊𝜕𝑎^,𝑝_𝜙=𝜕𝑊𝜕𝜙,𝑝^𝜔=𝜕𝑊𝜕𝜔^𝑝_𝐶1,=−𝜕𝑊𝜕𝐶₁,𝑝_𝐶2=−𝜕𝑊𝜕𝐶₂,𝑝_E=−𝜕𝑊𝜕E,(4.15) through Hamilton-Jacobi theory.

Let us consider the problem of determining a canonical transformation defined by a generating function 𝑊=𝑊(𝑎^,𝜙,𝜔^,𝐶₁,𝐶₂,E) such that 𝑎^2𝜇4𝑎^2𝜕𝜕𝑊𝑎^²+52𝜕𝜕𝑊𝜙²+52csc²𝜙−2𝜕𝑊𝜕𝜔^²=E.(4.16)𝑊 and 𝑊𝑎^,𝜙,𝜔^,𝐶₁,𝐶₂,E=𝑊¹𝑎^,𝐶₁,𝐶₂,E+𝑊²𝜙,𝐶¹,𝐶₂,E+𝑊³𝜔^,𝐶₁,𝐶₂,E.(4.17) become the new generalized coordinates, 𝜕𝑊₃𝜕𝜔^𝜕𝑊₂=𝐶₁,𝜕𝜙²+𝜕𝑊³𝜕𝜔^²𝑎_csc²𝜙=𝐶²²,2𝜇4𝑎^2𝜕𝜕𝑊𝑎^1²+52𝜕𝜕𝑊𝜙²²+52csc²𝜙−2𝜕𝑊³𝜕𝜔^²=E.(4.18)Since the new Hamiltonian function 𝑊₁=−5𝐶²24𝜇E5𝐶²𝑎^−1^1/2−tan⁻¹4𝜇E5𝐶²𝑎^𝑊₂−1^1/2,𝑊=𝐶²²−𝐶²¹csc²𝜙𝑑𝜙,₃=𝐶₁𝜔^,(4.19) 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 5𝐶²²−4𝐶²¹=5𝐶² where 𝑊₂ is the generating function. The corresponding Hamilton-Jacobi equation is then given by 𝐹^=EFollowing the separation of variables technique, it is assumed that the generating function 𝑑𝑝_𝐶1𝑑𝑑𝑝𝑡_𝐶2=0,𝑑𝑝𝑑_E𝑡=0,𝑑𝑡=−1,(4.20) is equal to a sum of functions, each of which depends on a single old variable, that is, 𝑝_𝐶1=𝛼₁,𝑝_𝐶2=𝛼₂,𝑝_E=𝛼₃−𝑡,(4.21)

Therefore, from (4.11) through (4.17), it follows that 𝛼_𝑖A complete solution of these equations is given by 𝑖=1,2,3 with 𝑊. We note that 𝑝^𝑎=12√2𝑎²𝑝_𝜙4𝜇E𝑎−5𝐶²𝑎²^1/2,𝑝=𝐶²²−𝐶²¹csc²𝜙,_𝜔𝛼₁=𝐶₁,=−𝜔−225𝐶𝐶¹tan⁻¹4𝜇E5𝐶²𝑎−1^1/2+𝐶¹csc²𝜙𝛼𝐶²²−𝐶²¹csc²𝜙𝑑𝜙,₂=−52𝐶𝐶²tan⁻¹4𝜇E5𝐶²𝑎−1^1/2+sin⁻¹𝐶²𝐶²²−𝐶²¹𝑡−𝛼₃cos𝜙,=−12√2E𝑎4𝜇E𝑎−5𝐶²𝑎²^1/2.(4.22) 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 Hamiltonian 𝐹^, that is, 𝑎^(𝑡)=𝑎^0𝑎1+4𝑎^0/𝜇(1/2)E𝑡²−𝑎^0𝑝^𝑎₀𝑡,(4.23)_𝜓s=in²1𝑘₀=𝑎^0sin²√2𝜓+𝑘₀,(4.24)5𝜏−𝜏⁰cos𝜙=cos𝑘1₁+4cos²𝑘₁,(4.25)𝜔cos𝜏,(4.26)_=𝑘₂+tan⁻¹tan𝜏csc𝑘¹−45𝜏sin𝑘¹𝑝,(4.27)_𝑎2=𝑎^0𝑎^³𝑝^𝑎2₀+18𝑝^𝜔205csc²𝑘₁−4𝑎^0𝑎^3−1𝑝𝑎^2,(4.28)_2𝜙=𝑝^𝜔2₀𝑝csc²𝑘₁−csc²𝜙,(4.29)_𝜔=𝑝^𝜔₀,(4.30) whose solution is very simple: 𝑘₀ where 𝑘₁,, 𝑘₂, are arbitrary constants of integration.

Introducing the generating function csc²𝑘₀=8𝑎^0𝑝^𝑎₀²+𝑝^𝜔2₀5csc²𝑘₁−4csc₂𝑝^𝜔2₀5csc²𝑘₁−4,𝑘₁=𝑝^2𝜙0+𝑝^𝜔2₀csc²𝜙₀𝑝^𝜔2₀𝑘₂,=𝜔^0+45𝜏⁰sin𝑘₁−tan⁻¹tan𝜏⁰csc𝑘₁.(4.31), 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 𝐶

From equations above, one finds the solution of the canonical system governed by the Hamiltonian 𝐶₁ for a given set of initial conditions: 𝐶₂,E𝐶²=15𝑝^𝜔205csc²𝐶𝑘₁−4₁,𝐶₂₂=𝑝^𝜔₀,4𝜇=𝑝^2𝜙₀+𝑝^𝜔2₀csc²𝜙₀,E=𝑎^08𝑎^0𝑝^𝑎₀²+𝑝^𝜔2₀5csc²𝑘₁−4.(4.32)𝑎^(0)=𝑎^0𝑒^(0)=sin𝜙₀,𝜔^(0)=𝜔^0𝜏₀cos𝜙₀=cos𝑘₁cos𝜏₀ with the auxiliary constants 𝑎𝑡=𝑎^𝑡+𝑎^5𝜇³8𝑒^sin𝐸^𝑎^2𝑝^𝑎+41−𝑒^2sin𝐸^𝑎^𝑝^𝑒−41−𝑒^2^1/2𝑒^cos𝐸^𝑎^𝑝^𝜔^𝐸′_𝐸′0𝑒𝑡=𝑒^𝑡+𝑎^5𝜇³41−𝑒^2sin𝐸^𝑎^𝑝^𝑎+1−𝑒^2−54𝑒^sin𝐸^+34sin2𝐸^−112𝑒^sin3𝐸^𝑝^𝑒+1−𝑒^2^1/2𝑒^54𝑒^cos𝐸^+14𝑒^2−3cos2𝐸^+112𝑒^𝜔𝑡=𝜔_cos3𝐸^𝑝^𝜔^𝐸′_𝐸′0,𝑡+𝑎^5𝜇³−41−𝑒^2^1/2𝑒^cos𝐸^𝑎^𝑝^𝑎+1−𝑒^2^1/2𝑒^54𝑒^cos𝐸^+14𝑒^2−3cos2𝐸^+112𝑒^cos3𝐸^𝑝^𝑒+1𝑒^245−𝑒^2𝑒^sin𝐸^−1232−𝑒^2sin2𝐸^+112𝑒^sin3𝐸^𝑝^𝜔^𝐸′_𝐸′0,(4.33), 𝑎^,𝑒^,…,𝑝^𝜔 and 𝑆 defined as functions of the initial value of the adjoint variables (conjugate momenta) by 𝐽The constants 𝐽=E𝑡−𝑡⁰+Δ𝑆,(4.34), Δ𝑆=𝑆(𝐸^)−𝑆(𝐸^0), cos²𝐸^ and cos2𝐸^ in (4.22) can also be written as functions of the initial value of the adjoint variables as follows: 𝐸^ The initial conditions for the state variables (generalized coordinates) are given by 𝑀^, 𝑀^(𝑡)=𝑀^0+^𝑡𝑡₀𝜇𝑎^3−5+2𝑒^22𝑎^5𝜇³√1−𝑒^2𝑒^2𝑝^𝜔𝑑𝑡,(4.35) and 𝑀^0=𝑀^(𝑡₀), and 𝑀^ is obtained from 𝐹₀.

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: Δ𝐹^ with 𝑝^𝑀 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 𝐹^=𝑛^−5+2𝑒^22𝑎^5𝜇³√1−𝑒^2𝑒^2𝑝^𝜔𝑝^𝑀+⋯,(4.36). The expression of 𝑝^𝜔₀=0 can be put in a compact form as follows: 𝑎^sin²𝑘₀=𝑎^0sin²(25𝜙−𝜙⁰+𝑘⁰𝑝),(4.37)_𝑎2=𝑎^0𝑎^³𝑝^𝑎2₀+58𝑝^2𝜙𝑎𝑎^0_3−1𝑝𝑎^2,(4.38)_𝜙=𝑝_𝜙0=𝐶,(4.39)with csc²𝑘₀=(8(𝑎^0𝑝^𝑎₀)²+5𝑝^2𝜙)/5𝑝^2𝜙.

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 4𝜇E=𝑎^0(8(𝑎^0𝑝^𝑎₀)²+5𝑝^2𝜙) appear instead of terms in 𝐻^∗; the other terms are equivalent. The eccentric anomaly 10⁻¹¹ is computed from Kepler's equation with the mean anomaly 10⁻¹² given by 10⁻⁴where 𝑎₀. The differential equation for the mean anomaly 𝑒₀ is derived from the undisturbed Hamiltonian 𝜔₀ and an additional term 𝑝_𝑎0 factored by 2.9326×10⁻⁴, 1.4663×10⁻⁴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. Since 𝑝_𝑒0, some equations must be rewritten: 2.9625×10⁻⁵1.4812×10⁻⁵𝑝_𝜔0 with 𝑡_𝑓−𝑡₀ and 𝑡_𝑓−𝑡₀.

Figures 3, 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 𝑎_𝑓, 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 𝑒_𝑓, and absolute error tolerance of 𝜔_𝑓, 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 −0.0094 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.

Table 2 

Final state variables.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Comparison between numerical, analytical, and average solutions for a transfer between coaxial orbits—set 1 of initial conditions and 𝜔.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Comparison between numerical, analytical, and average solutions for a transfer between coaxial orbits—set 2 of initial conditions and 𝜙.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

Comparison between numerical, analytical, and average solutions for a transfer between coaxial orbits—consumption variable.

5. Analysis of Long Duration Transfers

In this section, a complete analysis of long duration transfers described by the Hamiltonian 𝛼=𝛼(𝜙,𝜙₀,𝑘₀,𝑘₁) 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], 𝜔^=𝜔^(𝜙,𝜙₀,𝜔₀,𝑘₁), and (4.23) and (4.24) can be put in the form (𝑎^0,𝜙₀,𝜔^0)where 𝜕𝛼𝜕𝑘₀𝜕𝜔^𝜕𝑘₁−𝜕𝛼𝜕𝑘₁𝜕𝜔^𝜕𝑘₀=0.(5.5). By eliminating 𝜔^ from these equations, one finds 𝑘₀On the other hand, 𝜕𝜔^/𝜕𝑘₀=0, thus 𝜕𝛼𝜕𝑘₀=−2𝛼sin𝑘₀𝑣𝑢₀≠0(5.6)Note that 𝜕𝜔^𝜕𝑘₁=0.(5.7) in equations above.

Now, consider (4.24)–(4.27) which define a two-parameter family of extremals in the phase space 𝜔^=𝜔^(𝜙,𝜙₀,𝜔₀,𝑘₁) for a given initial phase point 𝜔^=𝜔^0−45sin𝑘¹cos⁻¹ccooss𝑘𝜙₁−cos⁻¹ccooss𝜙𝑘⁰₁+cos⁻¹ccoott𝑘𝜙₁−cos⁻¹ccoott𝜙𝑘⁰₁.(5.8) corresponding to an initial orbit. By eliminating the auxiliary variables 𝜕𝜔^/𝜕𝑘₁ and 𝜕𝜔^𝜕𝑘₁=−45cos𝑘¹𝜏−𝜏⁰−45tan²𝑘₁cos𝜙⁰sin𝜏₀−cos𝜙+sec₂sin𝜏𝑘₁[1√tan²𝜙₀−tan²𝑘₁−1√tan²𝜙−tan²𝑘₁].(5.9), 𝜏 and 1−45sin²𝑘₁sin𝜏−𝜏⁰−45cos²𝑘₁𝜏−𝜏⁰sin𝜏sin𝜏⁰=0.(5.10) can be written as explicit functions of 𝐹^, that is, 𝑒₀=0.4 and 𝑎₀=1. The conjugate points to the phase point 𝑘₀ are determined through the roots of the equation [12] 𝑢/𝑣₀Since 𝑒₀=0.0 does not depend on 𝑒₀=0.4, 𝑒₀=√3/2. On the other hand, from (4.24) and (5.2), one finds 𝑘₀=18^∘Thus, the problem of determining the conjugate points reduces to the analysis of the roots of the following equation: 𝑎₀=1.0

From (4.26) and (4.27), one finds the explicit form of 𝑒₀=0.5,, that is, 𝜔₀=0^∘The partial derivative 𝑘₁=18^∘ is given by 𝑘₀ 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]: 𝑘₀=18^∘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 (𝑒,𝜔) are globally minimizing.

Some numerical results are presented through Figures 6–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 𝑢/𝑣₀. In both cases, the semimajor axis is 𝑘₁=18^∘, such that the results are presented in canonical units. The extremals are plotted for different values of the constant [−180^∘,180^∘] and the isocost curves for different values of the new consumption variable 𝜔. Only transfers between direct orbits are considered.

Figure 6 

Field of extremals and isocost curves for transfers between coaxial orbits for 𝑂0∶(𝑎0,𝑒0,𝜔0).

Figure 7 

Field of extremals and isocost curves for transfers between coaxial orbits for 𝑂𝑓∶(𝑎𝑓,𝑒𝑓,𝜔𝑓).

(a)

(b)

Figure 8 

Field of extremals for transfers between noncoaxial orbits for 𝑘1 with 𝜙0,𝜙𝑓,𝜏0,.

Figure 9 

Isocost curves for 𝜏𝑓, 𝜔 and 𝜔≠𝜔𝑓, with 𝑘1.

Figures 8 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 |𝜔^−𝜔^𝑓|<𝛿 and 𝛿. In the plots of Figure 8, 𝜓_𝑓 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 (𝑢_𝑓/𝑣₀) -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).

Figure 9 shows a sample of isocost curves for different values of the new consumption variable 𝑘₀ and a fixed value of tan𝑘₀=sin(√2𝜓_𝑓)/(√𝛼_𝑓−cos(√2𝜓_𝑓)). Note that these curves lie on an isocost surface with complex geometry. For simplicity, the symmetric interval 𝛼_𝑓=𝑎^𝑓/𝑎^0. 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 orbit 𝑝^𝜔₀, to a final orbit 𝑝^𝜙₀ is presented. The steps of this algorithm can be described as follows.(1) Guess a starting value of 𝑝^𝑎₀=𝑣²⁰2𝑎^0𝑡_𝑓𝑢𝑣^𝑓₀cos𝑘⁰𝑝_𝜔0,=2√2𝑎^0𝑝^𝑎₀𝑝cot𝑘₀√5csc²𝑘₁−4_2𝜙0,=𝑝^𝜔2₀csc²𝑘₁−csc²𝜙₀.(5.11),(2) Determine 𝑘₁ and 𝑘_1𝑛+1=𝑘_1𝑛−𝜔^𝑘^1𝑛−𝜔^𝑓𝜕𝜔^/𝜕𝑘₁_{𝑘1=𝑘1𝑛},(5.12) through (4.9) and (4.26).(3) Compute (𝜕𝜔^/𝜕𝑘₁)_{𝑘1=𝑘1𝑛} through (4.27).(4) If 𝑘₁=0, adjust the value of 𝜓 and repeat steps (2) and (3) until (1/√5)(𝜙−𝜙₀), where 𝑝^𝑎₀ is a prescribed small quantity.(5) Compute 𝑝_𝜙0 and 𝑝^2𝜙₀=8tan²𝑘₀𝑎^0𝑝^𝑎₀²5.(5.13) through (4.25) and (5.3).(6) Compute 𝜸^∗=1𝑛𝑎√1−𝑒²2𝑎𝑒sin𝑓𝑝^𝑎+1−𝑒²sin𝑓𝑝^𝑒−1−𝑒²cos𝑓𝑒𝑝^𝜔𝐞^𝐫+2𝑎1−𝑒²𝑎𝑟𝑝^𝑎+1−𝑒²cos𝑓+cos𝐸𝑝^𝑒+1−𝑒²sin𝑓𝑒1+1+𝑒1cos𝑓𝑝^𝜔𝐞^𝐬.(5.14) using equation 𝑂₀∶(𝑎₀,𝑒₀,𝜔₀), where 𝑂_𝑓∶(𝑎_𝑓,𝑒_𝑓,𝜔_𝑓)(7) Compute successively 𝑡_𝑓, 𝑦_𝑖𝑡=𝑔_𝑖𝑡,𝑝^𝑎,𝑝^𝑒,𝑝^𝜔,𝑖=1,2,3,(6.1) and 𝑦₁(𝑡)=𝑎(𝑡) through the following equations: 𝑦₂(𝑡)=𝑒(𝑡), The adjust of the value of 𝑦₃(𝑡)=𝜔(𝑡) can be obtained through the well-known Newton-Raphson method, that is, 𝑝^𝑎 with the partial derivative 𝑝^𝑒, computed from (5.9).

In the case of coaxial orbits, 𝑝^𝜔 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 𝑒^. The initial value of the adjoint variables is given as follows: 𝜔^, is computed through step (7) and 𝐸^ is given by 𝑔_𝑖

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, 𝑖=1,2,3 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 𝑝^𝑎 to a final orbit 𝑝^𝑒, at the prescribed final time 𝑝^𝜔.

For a given time, (4.33) can be represented as follows: 𝑦₁(𝑡_𝑓)=𝑎_𝑓where 𝑦₂(𝑡_𝑓)=𝑒_𝑓,, 𝑦₃(𝑡_𝑓)=𝜔_𝑓 and 𝑂₀. Note that 𝑎₀=1.0, 𝑒₀=0.20 and 𝜔₀=0^∘ appear explicitly in the short periodic terms, but they also appear implicitly through 𝑎_𝑓=2.0, 𝑒_𝑓=0.25, 𝜔_𝑓=0^∘ and 𝑡_𝑓−𝑡₀=500. Thus, functions 𝑡_𝑓−𝑡₀=1000, 𝑡_𝑓−𝑡₀=500, are nonlinear in these variables.

So, the two-point boundary value problem can be stated as follows: find 𝑡_𝑓−𝑡₀=1000, 𝐻^∗ and 𝐹₁ such that 𝑆₁, 𝑂 and 𝑎. 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.

Figure 10 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 𝑒, defined by the set of initial conditions 𝜔, 𝑎, 𝑒, and the arrival orbit defined by the set of terminal conditions 𝜔, 𝑎^, 𝑒^, and two transfer durations 𝜔^ and Δ𝑥=𝐴𝑝₀,(7.1). All variables are represented in canonical units.

(a)

(b)

Figure 10 

Complete and average solutions of the two-point boundary value problem for Δ𝑥=[Δ𝛼Δ𝑒Δ𝜔]𝑇 and 𝛼=𝑎/𝑎.

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 Hamiltonian 𝑝_𝛼=𝑎𝑝_𝑎, defined through (3.21) and (3.22), can be straightforwardly obtained by linearizing the functions 𝑝₀ and 3×1 about a reference elliptic orbit 𝐴 defined by a nominal set of orbital elements 3×3, 𝐴 and 𝐴=⎡⎢⎢⎢⎣𝑎^𝛼𝛼𝑎_𝛼𝑒𝑎_𝛼𝜔𝑎_𝑒𝛼𝑎_𝑒𝑒𝑎_𝑒𝜔𝑎_𝜔𝛼𝑎_𝜔𝑒𝑎_𝜔𝜔⎤⎥⎥⎥⎦,(7.2). This solution is very similar to the complete solution given by (4.33), with 𝑎_𝛼𝛼=4𝑎⁵𝜇³𝑎_𝛼𝑒𝐸+𝑒sin𝐸|||^𝐸𝑓_𝐸0,=𝑎_𝑒𝛼=4𝑎⁵𝜇³1−𝑎_𝛼𝜔𝑒²sin𝐸|||^𝐸𝑓_𝐸0,=𝑎_𝜔𝛼=−4𝑎⁵𝜇³1−𝑒²^1/2𝑒cos𝑎_𝑒𝑒𝐸|||^𝐸𝑓_𝐸0,=𝑎⁵𝜇³1−𝑒²52𝑀−54𝑒sin𝐸+34sin2𝐸−112𝑎_𝑒𝜔𝑒sin3𝐸|||^𝐸𝑓_𝐸0,=𝑎_𝜔𝑒=𝑎⁵𝜇³1−𝑒²^1/2𝑒45𝑒cos𝐸+14𝑒²−3cos2𝐸+112𝑎_𝜔𝜔𝑒cos3𝐸|||^𝐸𝑓_𝐸0,=𝑎⁵𝜇³1𝑒²52−2𝑒²𝑀+𝑒54−𝑒²sin𝐸−1223−𝑒²sin2𝐸+112𝑒sin3𝐸|||^𝐸𝑓_𝐸0,(7.3), 𝐸 and 𝑀=𝑛(𝑡−𝜏) replacing 𝑛=𝜇/𝑎³, 𝜏, 𝑂₀∶(𝑎₀,𝑒₀,𝜔₀), and can be put in a suitable matrix form 𝑂_𝑓∶(𝑎_𝑓,𝑒_𝑓,𝜔_𝑓)where 𝑡_𝑓 denotes the imposed changes on orbital elements (state variable), 𝑝₀, 𝑡_𝑓−𝑡₀=5, 𝑎 is the 𝑒 matrix of initial values of the adjoint variable and 𝜔 is 𝐹₁ symmetric matrix. The adjoint variables are constant, and the matrix 𝑆₁ is given by 𝑂 where 𝐻(𝑥,𝑦,𝜀) where 𝜀 is calculated by solving Kepler’s equation for 𝐻(𝑥,𝑦,𝜀)=𝐻₀(𝑥,𝑦)+_𝑘=1𝜀^𝑘𝐻_𝑘(𝑥,𝑦),(A.1), with 𝐻₀(𝑥,𝑦). 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 𝑦 to a final orbit 𝑛 at the prescribed final time (𝑥,𝑦)→(𝑥^,𝑦^). In this case, the vector 𝑆(𝑥^,𝑦^,𝜀) 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.

Figure 11 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 𝜀 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11 

Linearized and nonlinear solutions for transfers between close orbits.

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 function 𝑛 and the generating function 𝐹(𝑥^,𝑦^,𝜀) about a reference elliptic orbit 𝜀. 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.

Appendix

Hori method is based on Lie theorem [2].

Let 𝐹𝑥^,𝑦^=𝐹⁰𝑥^,𝑦^+_𝑘=1𝜀^𝑘𝐹_𝑘𝑥^,𝑦^.(A.3) be a Hamiltonian written as a power series of a small parameter 𝑥_𝑖=𝑥^𝑖+^∞_𝑛=1𝜀^𝑛𝑛!𝐷^𝑆𝑛−1𝜕𝑆𝜕𝑦^𝑖𝑦_𝑖,=𝑦^𝑖+^∞_𝑛=1𝜀^𝑛𝑛!𝐷^𝑆𝑛−1𝜕𝑆𝜕𝑥^𝑖,𝑖=1,…,𝑛,(A.4): 𝐷^0𝑆𝑓=𝑓where the undisturbed Hamiltonian 𝐷^1𝑆𝑓={𝑓,𝑆} describes an integrable system. Here, 𝐷^𝑛𝑆𝑓=𝐷^𝑆𝑛−1(𝐷_𝑆𝑓) and 𝑛≥2 are 𝑓(𝑥^,𝑦^)-vectors of generalized coordinates and conjugate momenta, respectively.

Consider an infinitesimal canonical transformation 𝐹(𝑥^,𝑦^,𝜀) defined by a generating function 𝑆(𝑥^,𝑦^,𝜀), also developed as a power series of the small parameter 𝜀, 𝐻₀𝑥^,𝑦^=𝐹⁰𝑥^,𝑦^,(A.5)This transformation is such that the new dynamical system has some advantages for the solution. Here, prime denotes the new variables, and 𝐻⁰,𝑆₁+𝐻¹=𝐹₁,(A.6) and 𝐻⁰,𝑆₂+21𝐻¹+𝐹₁,𝑆₁+𝐻²=𝐹₂,….(A.7) are also (𝑥^,𝑦^)-vectors. The new Hamiltonian 𝑚, resulting from this canonical transformation, is also developed as a power series of the small parameter 𝐻⁰,𝑆_𝑚+Θ^𝑚=𝐹_𝑚,(A.8), Θ_𝑚

The transformation equations are given explicitly by 𝐻₀ where 𝐻_𝑚, 𝑆_𝑘, 𝐹_𝑘,, 𝐻_𝑘, 𝑘=1,…,𝑚−1 denotes an arbitrary function and braces stand for Poisson brackets.

According to the algorithm of Hori method, the new Hamiltonian 𝐹_𝑚 and the generating function 𝑆_𝑚 are obtained, at each order of the small parameter 𝑑𝑥^𝑖𝑑𝑡=𝜕𝜕𝐹𝑦^0𝑖,𝑑𝑦^𝑖𝑑𝑡=−𝜕𝜕𝐹𝑥^0𝑖,𝑖=1,…,𝑛,(A.9), from the equations:(i)order 0: 𝜕𝑆_𝑚𝜕𝑡=Θ^𝑚−𝐹_𝑚,(A.10)(ii)order 1: Θ_𝑚(iii)order 2: 𝑐₁,…,𝑐_2𝑛 All functions in above equations are written in terms of the new set of canonical variables 𝑆_𝑚.

The 𝐹_𝑚th order equation of the algorithm can be put in the form 𝑆where 𝑡 is a function obtained from the preceding orders, involving 𝐹, 𝐹_𝑚, Θ_𝑚, 𝑆_𝑚 and 𝐹_𝑚, 𝑆𝐹_𝑚𝑚=Θ_𝑚,=Θ_𝑚−Θ_𝑚𝑑𝑡,(A.11).

The determination of the functions ⟨⋅⟩ and is based on the general solution of the undisturbed system and it is performed through an averaging principle. Following the integration theory proposed in [18], (A.8) can be put in the form with written in terms of the arbitrary constants of integration, , of the general solution of the canonical system (A.9). and are unknown functions.

Following Hori [2], suppose that the canonical transformation generated by the function is such that the time is eliminated from the new Hamiltonian . This is accomplished taking as the mean value of the function . Therefore, functions and are given through the equations where stands for the mean value of the function. It should be noted that the averaging process and the integration in (A.11) are performed considering the explicit dependence on time of the functions.

Acknowledgment

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