Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2008 (2008), Article ID 525930, 30 pages
http://dx.doi.org/10.1155/2008/525930
Research Article

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

Departamento de Matemática, Instituto Tecnológico de Aeronáutica, 12228-900 São José dos Campos, SP, Brazil

Received 7 March 2008; Accepted 15 May 2008

Academic Editor: Alexander P. Seyranian

Copyright © 2008 Sandro da Silva Fernandes and Francisco das Chagas Carvalho. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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/𝑔0, is between 104 and 102. 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 the space vehicle, 𝐽=𝑃max1𝑚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 𝑚𝑓.

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,𝐯0,0) at time 𝑡0 to the final state (𝐫𝑓,𝐯𝑓,𝐽𝑓) at time 𝑡𝑓, such that the final consumption variable 𝐽𝑓 is a minimum (Figure 1).

525930.fig.001
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 𝑡𝑓𝑡0 and the position of the vehicle in the initial orbit are specified. The state equations are 𝑑𝐫𝑑𝑡=𝐯,𝑑𝐯𝜇𝑑𝑡=𝑟3𝐫+𝜸,𝑑𝐽=1𝑑𝑡2𝛾2,(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), 𝐻=𝐩𝑟𝐯+𝐩𝑣𝜇𝑟3+1𝐫+𝜸2𝑝𝐽𝛾2,(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 𝐻, 𝐻=𝐩𝑟𝐯𝐩𝑣𝜇𝑟3𝑝𝐫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, 𝐻=𝐩𝑟𝐯𝐩𝑣𝜇𝑟3𝑝𝐫+2𝑣2.(2.7)

This Hamiltonian can be decomposed into two parts: 𝐻0=𝐩𝑟𝐯𝐩𝑣𝜇𝑟3𝐻𝐫,𝛾=𝑝2𝑣2,(2.8)where 𝐻0 is the undisturbed Hamiltonian and 𝐻𝛾 is the disturbing function concerning the optimal thrust acceleration. The undisturbed Hamiltonian 𝐻0 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 𝐻0, 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𝑣𝐞𝐫𝐞𝐫,𝑑𝐯𝜇𝑑𝑡=𝑟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𝑒2𝐞1+𝑒cos𝑓𝐫,(3.2)𝐯=𝜇𝑎1𝑒2𝑒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𝑟𝑣2𝑒/𝜇,(3.5)2=12𝜇𝑎,(3.6)cos𝜔=𝐢𝐞𝑒1,(3.7)cos𝐸=𝑒𝑟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𝐞=𝜇𝑣2𝜇𝑟𝐯,𝐫𝐫𝐯𝐡=𝐫×𝐯.(3.9) and the angular momentum vector ×, shown in Figure 2, are given, as function of the Cartesian elements 𝑓 and 𝐸, by the following equations: 𝑓tan2=1+𝑒𝐸1𝑒tan2.(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 𝛿𝐯

525930.fig.002
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𝑎2𝜇𝐯𝛿𝐯,𝛿𝑒=𝐡𝛿𝐡+𝜇𝑎𝑒22𝜇𝑎2𝑒𝛿𝑎,𝛿𝜔=𝛿𝑒𝑒cot𝜔𝐢𝛿𝐞,1𝑒sin𝜔𝛿𝐸=𝑒sin𝐸𝐫𝛿𝐫𝑟𝑟𝑎𝑎2,𝛿𝑎+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𝛿𝜉+2𝑒sin𝑓𝑛1𝑒22𝛿𝑢+1𝑒2𝑛𝑎𝑟𝑎𝛿𝑣,𝛿𝑒=1𝑒2𝑟2cos𝐸𝛿𝜉+sin𝑓𝑎𝛿𝜂+1𝑒2𝑛𝑎sin𝑓𝛿𝑢+1𝑒2𝑛𝑎cos𝐸+cos𝑓𝛿𝑣,𝛿𝜔=sin𝑓𝑒𝑟𝛿𝜉𝑒+cos𝑓𝑎𝑒1𝑒2𝛿𝜂1𝑒2𝑛𝑎𝑒cos𝑓𝛿𝑢+1𝑒21𝑛𝑎𝑒sin𝑓1+1+𝑒cos𝑓𝛿𝑣,𝛿𝑀=1𝑒3cos𝐸𝑒𝑟1𝑒2sin𝑓𝛿𝜉+1𝑒2+𝑎𝑒cos𝑓𝛿𝜂1𝑒2𝑛𝑎𝑒cos𝑓2𝑒1+𝑒cos𝑓𝛿𝑢1𝑒21𝑛𝑎𝑒sin𝑓1+1+𝑒cos𝑓𝛿𝑣,(3.17). At this point, it is useful to replace the variation in the eccentric anomaly 𝑛=𝜇/𝑎3 by the variation in the mean anomaly 𝛿𝑎𝛿𝑒𝛿𝜔𝛿𝑀=Δ1𝛿𝜉𝛿𝜂𝛿𝑢𝛿𝑣,(3.18), obtained from the well-known Kepler equation Δ1This variation is given by 𝐩𝐫Thus, 𝐩𝐯 where 𝐩𝐫𝐩𝐯=Δ1𝑇𝑝𝑎𝑝𝑒𝑝𝜔𝑝𝑀,(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 𝐩𝐫=𝑎𝑟22𝑎𝑝𝑎+1𝑒2𝑝cos𝐸𝑒+𝑟𝑎sin𝑓𝑒𝑝𝜔1𝑒3cos𝐸1𝑒2𝑝𝑀𝐞𝐫+sin𝑓𝑎𝑝𝑒𝑒+cos𝑓𝑎𝑒1𝑒2𝑝𝜔+1𝑒2cos𝑓𝑝𝑎𝑒𝑀𝐞𝐬,𝐩𝐯=1𝑛𝑎1𝑒22𝑎𝑒sin𝑓𝑝𝑎+1𝑒2𝑝sin𝑓𝑒1𝑒2cos𝑓𝑒𝑝𝜔+1𝑒23/2𝑒cos𝑓2𝑒𝑝1+𝑒cos𝑓𝑀𝐞𝐫+2𝑎1𝑒2𝑎𝑟𝑝𝑎+1𝑒2𝑝cos𝑓+cos𝐸𝑒+1𝑒2sin𝑓𝑒11+1+𝑒cos𝑓(𝑝𝜔1𝑒2𝑝𝑀)𝐞𝐬.(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 𝐻𝐻0=𝑛𝑝𝑀,𝐻(3.21)𝛾=12𝑛2𝑎21𝑒2×12(1cos2𝑓)2𝑎𝑒𝑝𝑎+1𝑒2𝑝𝑒2+21𝑒2sin2𝑓𝑎𝑝𝑎𝑝𝜔1𝑒2𝑝2𝑒𝑒𝑝𝜔+41𝑒23/2sin𝑓2𝑒1+𝑒cos𝑓+cos𝑓𝑎𝑝𝑎𝑝𝑀+1𝑒2𝑝2𝑒𝑒𝑝𝑀+1𝑒222𝑒2(1+cos2𝑓)𝑝2𝜔21𝑒25/2𝑒22𝑒1+𝑒cos𝑓+cos𝑓cos𝑓𝑝𝜔𝑝𝑀+1𝑒23𝑒22𝑒1+𝑒cos𝑓+cos𝑓2𝑝2𝑀+4𝑎21𝑒22𝑎𝑟2𝑝2𝑎+4𝑎1𝑒22𝑎𝑟(cos𝐸+cos𝑓)𝑝𝑎𝑝𝑒+1𝑒22(cos𝐸+cos𝑓)2𝑝2𝑒+4𝑎1𝑒22𝑒𝑎𝑟1sin𝑓1+1+𝑒cos𝑓[𝑝𝑎𝑝𝜔1𝑒2𝑝𝑎𝑝𝑀]+21𝑒22𝑒1(cos𝐸+cos𝑓)1+1+𝑒cos𝑓sin𝑓[𝑝𝑒𝑝𝜔1𝑒2𝑝𝑒𝑝𝑀]+1𝑒2𝑒11+1+𝑒cos𝑓sin𝑓[𝑝𝜔1𝑒2𝑝𝑀]2.(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 𝐻0 is invariant with respect to this canonical transformation. Thus, 𝐻𝛾

The new Hamiltonian function 𝐻0, 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 𝐹0=𝑛𝑝𝑀.(4.2) and 𝐹0 have different order of magnitude: 𝑑𝑎𝑑𝑡=0,𝑑𝑝𝑎=3𝑑𝑡2𝑛𝑎𝑝𝑀,𝑑𝑒𝑑𝑡=0,𝑑𝑝𝑒𝑑𝑡=0,𝑑𝜔𝑑𝑡=0,𝑑𝑝𝜔𝑑𝑡=0,𝑑𝑀𝑑𝑡=𝑛,𝑑𝑝𝑀𝑑𝑡=0,(4.3) has zero order and 𝑎=𝑎0,𝑒=𝑒0,𝜔=𝜔0,𝑀=𝑀0+𝑛𝑡𝑡0,𝑝𝑎=𝑝𝑎0+32𝑛𝑡𝑡0𝑎𝑝𝑀,𝑝𝑒=𝑝𝑒0,𝑝𝜔=𝑝𝜔0,𝑝𝑀=𝑝𝑀0.(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, 𝜕𝑆1𝜕𝑡=𝐻𝛾𝐹1.(4.5) Following the algorithm described in the appendix, at order 0 one gets 𝐻𝛾

Now, consider the undisturbed canonical system described by 𝐹1 : 𝐻𝛾 general solution of which is given by 𝐹1=𝑎2𝜇4𝑎2𝑝2𝑎+521𝑒2𝑝2𝑒+54𝑒22𝑒2𝑝2𝜔,𝑆(4.6)1=12𝑎5𝜇38𝑒sin𝐸𝑎2𝑝𝑎2+81𝑒2sin𝐸𝑎𝑝𝑎𝑝𝑒81𝑒21/2𝑒cos𝐸𝑎𝑝𝑎𝑝𝜔+1𝑒254𝑒sin𝐸+34sin2𝐸1𝑒12sin3𝐸𝑝𝑒2+21𝑒21/2𝑒54𝑒cos𝐸+14𝑒23cos2𝐸+1𝑒12cos3𝐸𝑝𝑒𝑝𝜔+1𝑒254𝑒2𝑒sin𝐸1232𝑒2sin2𝐸+1𝑒12sin3𝐸𝑝𝜔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 𝐹=𝐹1 and 𝑆=𝑆1 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+52csc2𝑝𝜙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 𝑝𝜔=𝐶1𝑝,(4.11)2𝜙+𝑝𝜔2csc2𝜙=𝐶22𝑎,(4.12)2𝜇4𝑎2𝑝𝑎2+52𝑝𝜙2+52csc2𝑝𝜙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), 𝐶1,𝐶2,E(Constant 𝑎,𝜙,𝜔,𝑝𝑎,𝑝𝜙,𝑝𝜔𝑊𝐶1,𝐶2,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=𝜕𝑊𝜕𝐶1,𝑝𝐶2=𝜕𝑊𝜕𝐶2,𝑝E=𝜕𝑊𝜕E,(4.15) through Hamilton-Jacobi theory.

Let us consider the problem of determining a canonical transformation defined by a generating function 𝑊=𝑊(𝑎,𝜙,𝜔,𝐶1,𝐶2,E) such that 𝑎2𝜇4𝑎2𝜕𝑊𝜕𝑎2+52𝜕𝑊𝜕𝜙2+52csc2𝜙2𝜕𝑊𝜕𝜔2=E.(4.16)𝑊 and 𝑊𝑎,𝜙,𝜔,𝐶1,𝐶2,E=𝑊1𝑎,𝐶1,𝐶2,E+𝑊2𝜙,𝐶1,𝐶2,E+𝑊3𝜔,𝐶1,𝐶2,E.(4.17) become the new generalized coordinates, 𝜕𝑊3𝜕𝜔=𝐶1,𝜕𝑊2𝜕𝜙2+𝜕𝑊3𝜕𝜔2csc2𝜙=𝐶22,𝑎2𝜇4𝑎2𝜕𝑊1𝜕𝑎2+52𝜕𝑊2𝜕𝜙2+52csc2𝜙2𝜕𝑊3𝜕𝜔2=E.(4.18)Since the new Hamiltonian function 𝑊1=5𝐶224𝜇E5𝐶2𝑎11/2tan14𝜇E5𝐶2𝑎11/2,𝑊2=𝐶22𝐶21csc2𝑊𝜙𝑑𝜙,3=𝐶1𝜔,(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𝐶224𝐶21=5𝐶2 where 𝑊2 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𝑑𝑡=0,𝑑𝑝𝐶2𝑑𝑡=0,𝑑𝑝E𝑑𝑡=1,(4.20) is equal to a sum of functions, each of which depends on a single old variable, that is, 𝑝𝐶1=𝛼1,𝑝𝐶2=𝛼2,𝑝E=𝛼3𝑡,(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 𝑝𝑎=122𝑎24𝜇E𝑎5𝐶2𝑎21/2,𝑝𝜙=𝐶22𝐶21csc2𝑝𝜙,𝜔=𝐶1,𝛼1=𝜔225𝐶1𝐶tan14𝜇E5𝐶2𝑎11/2+𝐶1csc2𝜙𝐶22𝐶21csc2𝜙𝛼𝑑𝜙,2=52𝐶2𝐶tan14𝜇E5𝐶2𝑎11/2+sin1𝐶2𝐶22𝐶21,cos𝜙𝑡𝛼31=22E𝑎4𝜇E𝑎5𝐶2𝑎21/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, 𝑎𝑎(𝑡)=01+4𝑎0/𝜇(1/2)E𝑡2𝑎0𝑝𝑎0𝑡𝑎,(4.23)sin2𝑘0=𝑎0sin22𝜓+𝑘01,(4.24)𝜓=5𝜏𝜏01+4cos2𝑘1,(4.25)cos𝜙=cos𝑘1𝜔cos𝜏,(4.26)=𝑘2+tan1tan𝜏csc𝑘145𝜏sin𝑘1𝑝,(4.27)𝑎2=𝑎0𝑎3𝑝𝑎20+18𝑝𝜔205csc2𝑘1𝑎40𝑎31𝑎2𝑝,(4.28)2𝜙=𝑝𝜔20csc2𝑘1csc2𝜙𝑝,(4.29)𝜔=𝑝𝜔0,(4.30) whose solution is very simple: 𝑘0 where 𝑘1,, 𝑘2, are arbitrary constants of integration.

Introducing the generating function csc2𝑘0=8𝑎0𝑝𝑎02+𝑝𝜔205csc2𝑘14𝑝𝜔205csc2𝑘1,4csc2𝑘1=𝑝2𝜙0+𝑝𝜔20csc2𝜙0𝑝𝜔20,𝑘2=𝜔0+45𝜏0sin𝑘1tan1tan𝜏0csc𝑘1.(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 𝐶1 for a given set of initial conditions: 𝐶2,E𝐶2=15𝑝𝜔205csc2𝑘1,𝐶41=𝑝𝜔0,𝐶22=𝑝2𝜙0+𝑝𝜔20csc2𝜙0,4𝜇E=𝑎08𝑎0𝑝𝑎02+𝑝𝜔205csc2𝑘1.4(4.32)𝑎(0)=𝑎0𝑒(0)=sin𝜙0,𝜔(0)=𝜔0𝜏0cos𝜙0=cos𝑘1cos𝜏0 with the auxiliary constants 𝑎𝑡=𝑎𝑡+𝑎5𝜇38𝑒sin𝐸𝑎2𝑝𝑎+41𝑒2sin𝐸𝑎𝑝𝑒41𝑒21/2𝑒cos𝐸𝑎𝑝𝜔𝐸𝐸0𝑒𝑡=𝑒𝑡+𝑎5𝜇341𝑒2sin𝐸𝑎𝑝𝑎+1𝑒254𝑒sin𝐸+34sin2𝐸1𝑒12sin3𝐸𝑝𝑒+1𝑒21/2𝑒54𝑒cos𝐸+14𝑒23cos2𝐸+1𝑒12cos3𝐸𝑝𝜔𝐸𝐸0,𝜔𝑡=𝜔𝑡+𝑎5𝜇341𝑒21/2𝑒cos𝐸𝑎𝑝𝑎+1𝑒21/2𝑒54𝑒cos𝐸+14𝑒23cos2𝐸+1𝑒12cos3𝐸𝑝𝑒+1𝑒254𝑒2𝑒sin𝐸1232𝑒2sin2𝐸+1𝑒12sin3𝐸𝑝𝜔𝐸𝐸0,(4.33), 𝑎,𝑒,,𝑝𝜔 and 𝑆 defined as functions of the initial value of the adjoint variables (conjugate momenta) by 𝐽The constants 𝐽=E𝑡𝑡0+Δ𝑆,(4.34), Δ𝑆=𝑆(𝐸)𝑆(𝐸0), cos2𝐸 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+𝑡𝑡0𝜇𝑎35+2𝑒22𝑎5𝜇31𝑒2𝑒2𝑝𝜔𝑑𝑡,(4.35) and 𝑀0=𝑀(𝑡0), and 𝑀 is obtained from 𝐹0.

Equation (4.23) through (4.30) represents the solution of the canonical system concerning the problem of optimal long duration low-thrust limited-power transfers between arbitrary elliptic coplanar orbits in a Newtonian central field.

For maneuvers with arbitrary duration, the periodic terms must be included [9]. Following Hori method and applying the initial conditions, one gets a first-order analytical solution: Δ𝐹 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𝜇31𝑒2𝑒2𝑝𝜔𝑝𝑀+,(4.36). The expression of 𝑝𝜔0=0 can be put in a compact form as follows: 𝑎sin2𝑘0=𝑎0sin2(25𝜙𝜙0+𝑘0𝑝),(4.37)𝑎2=𝑎0𝑎3𝑝𝑎20+58𝑝2𝜙𝑎0𝑎31𝑎2𝑝,(4.38)𝜙=𝑝𝜙0=𝐶,(4.39)with csc2𝑘0=(8(𝑎0𝑝𝑎0)2+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𝑝𝑎0)2+5𝑝2𝜙) appear instead of terms in 𝐻; the other terms are equivalent. The eccentric anomaly 1011 is computed from Kepler's equation with the mean anomaly 1012 given by 104where 𝑎0. The differential equation for the mean anomaly 𝑒0 is derived from the undisturbed Hamiltonian 𝜔0 and an additional term 𝑝𝑎0 factored by 2.9326×104, 1.4663×104where 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×1051.4812×105𝑝𝜔0 with 𝑡𝑓𝑡0 and 𝑡𝑓𝑡0.

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.

tab1
Table 1: Set of initial conditions and transfer duration.
tab2
Table 2: Final state variables.
fig3
Figure 3: Comparison between numerical, analytical, and average solutions for a transfer between coaxial orbits—set 1 of initial conditions and 𝜔.
fig4
Figure 4: Comparison between numerical, analytical, and average solutions for a transfer between coaxial orbits—set 2 of initial conditions and 𝜙.
fig5
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 𝛼=𝛼(𝜙,𝜙0,𝑘0,𝑘1) 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], 𝜔=𝜔(𝜙,𝜙0,𝜔0,𝑘1), and (4.23) and (4.24) can be put in the form (𝑎0,𝜙0,𝜔0)where 𝜕𝛼𝜕𝑘0𝜕𝜔𝜕𝑘1𝜕𝛼𝜕𝑘1𝜕𝜔𝜕𝑘0=0.(5.5). By eliminating 𝜔 from these equations, one finds 𝑘0On the other hand, 𝜕𝜔/𝜕𝑘0=0, thus 𝜕𝛼𝜕𝑘0=2𝛼sin𝑘0𝑢𝑣00(5.6)Note that 𝜕𝜔𝜕𝑘1=0.(5.7) in equations above.

Now, consider (4.24)–(4.27) which define a two-parameter family of extremals in the phase space 𝜔=𝜔(𝜙,𝜙0,𝜔0,𝑘1) for a given initial phase point 𝜔=𝜔045sin𝑘1cos1cos𝜙cos𝑘1cos1cos𝜙0cos𝑘1+cos1cot𝜙cot𝑘1cos1cot𝜙0cot𝑘1.(5.8) corresponding to an initial orbit. By eliminating the auxiliary variables 𝜕𝜔/𝜕𝑘1 and 𝜕𝜔𝜕𝑘14=5cos𝑘1𝜏𝜏045tan2𝑘1cos𝜙0sin𝜏0cos𝜙sin𝜏+sec2𝑘1[1tan2𝜙0tan2𝑘11tan2𝜙tan2𝑘1].(5.9), 𝜏 and 415sin2𝑘1sin𝜏𝜏045cos2𝑘1𝜏𝜏0sin𝜏sin𝜏0=0.(5.10) can be written as explicit functions of 𝐹, that is, 𝑒0=0.4 and 𝑎0=1. The conjugate points to the phase point 𝑘0 are determined through the roots of the equation [12] 𝑢/𝑣0Since 𝑒0=0.0 does not depend on 𝑒0=0.4, 𝑒0=3/2. On the other hand, from (4.24) and (5.2), one finds 𝑘0=18Thus, the problem of determining the conjugate points reduces to the analysis of the roots of the following equation: 𝑎0=1.0

From (4.26) and (4.27), one finds the explicit form of 𝑒0=0.5,, that is, 𝜔0=0The partial derivative 𝑘1=18 is given by 𝑘0 The auxiliary variable 𝑘1 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]: 𝑘0=18As 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 69. 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 𝑢/𝑣0. In both cases, the semimajor axis is 𝑘1=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.

525930.fig.006
Figure 6: Field of extremals and isocost curves for transfers between coaxial orbits for 𝑂0(𝑎0,𝑒0,𝜔0).
525930.fig.007
Figure 7: Field of extremals and isocost curves for transfers between coaxial orbits for 𝑂𝑓(𝑎𝑓,𝑒𝑓,𝜔𝑓).
fig8
Figure 8: Field of extremals for transfers between noncoaxial orbits for 𝑘1 with 𝜙0,𝜙𝑓,𝜏0,.
525930.fig.009
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 (𝑢𝑓/𝑣0) -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 𝑘0 and a fixed value of tan𝑘0=sin(2𝜓𝑓)/(𝛼𝑓cos(2𝜓𝑓)). Note that these curves lie on an isocost surface with complex geometry. For simplicity, the symmetric interval 𝛼𝑓=𝑎𝑓/𝑎0. for 𝑝𝑎0 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 𝑝𝜔0, to a final orbit 𝑝𝜙0 is presented. The steps of this algorithm can be described as follows.(1) Guess a starting value of 𝑝𝑎0=𝑣202𝑎0𝑡𝑓𝑢𝑓𝑣0cos𝑘0,𝑝𝜔0=22𝑎0𝑝𝑎0cot𝑘05csc2𝑘1,𝑝42𝜙0=𝑝𝜔20csc2𝑘1csc2𝜙0.(5.11),(2) Determine 𝑘1 and 𝑘1𝑛+1=𝑘1𝑛𝜔𝑘1𝑛𝜔𝑓𝜕𝜔/𝜕𝑘1𝑘1=𝑘1𝑛,(5.12) through (4.9) and (4.26).(3) Compute (𝜕𝜔/𝜕𝑘1)𝑘1=𝑘1𝑛 through (4.27).(4) If 𝑘1=0, adjust the value of 𝜓 and repeat steps (2) and (3) until (1/5)(𝜙𝜙0), where 𝑝𝑎0 is a prescribed small quantity.(5) Compute 𝑝𝜙0 and 𝑝2𝜙0=8tan2𝑘0𝑎0𝑝𝑎025.(5.13) through (4.25) and (5.3).(6) Compute 𝜸=1𝑛𝑎1𝑒22𝑎𝑒sin𝑓𝑝𝑎+1𝑒2𝑝sin𝑓𝑒1𝑒2cos𝑓𝑒𝑝𝜔𝐞𝐫+2𝑎1𝑒2𝑎𝑟𝑝𝑎+1𝑒2𝑝cos𝑓+cos𝐸𝑒+1𝑒2sin𝑓𝑒11+𝑝1+𝑒cos𝑓𝜔𝐞𝐬.(5.14) using equation 𝑂0(𝑎0,𝑒0,𝜔0), where 𝑂𝑓(𝑎𝑓,𝑒𝑓,𝜔𝑓)(7) Compute successively 𝑡𝑓, 𝑦𝑖𝑡=𝑔𝑖𝑡,𝑝𝑎,𝑝𝑒,𝑝𝜔,𝑖=1,2,3,(6.1) and 𝑦1(𝑡)=𝑎(𝑡) through the following equations: 𝑦2(𝑡)=𝑒(𝑡), The adjust of the value of 𝑦3(𝑡)=𝜔(𝑡) 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: 𝑦1(𝑡𝑓)=𝑎𝑓where 𝑦2(𝑡𝑓)=𝑒𝑓,, 𝑦3(𝑡𝑓)=𝜔𝑓 and 𝑂0. Note that 𝑎0=1.0, 𝑒0=0.20 and 𝜔0=0 appear explicitly in the short periodic terms, but they also appear implicitly through 𝑎𝑓=2.0, 𝑒𝑓=0.25, 𝜔𝑓=0 and 𝑡𝑓𝑡0=500. Thus, functions 𝑡𝑓𝑡0=1000, 𝑡𝑓𝑡0=500, are nonlinear in these variables.

So, the two-point boundary value problem can be stated as follows: find 𝑡𝑓𝑡0=1000, 𝐻 and 𝐹1 such that 𝑆1, 𝑂 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 Δ𝑥=𝐴𝑝0,(7.1). All variables are represented in canonical units.

fig10
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 𝑝0 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𝑎5𝜇3𝐸+𝑒sin𝐸|||𝐸𝑓𝐸0,𝑎𝛼𝑒=𝑎𝑒𝛼=4𝑎5𝜇31𝑒2sin𝐸|||𝐸𝑓𝐸0,𝑎𝛼𝜔=𝑎𝜔𝛼=4𝑎5𝜇31𝑒21/2𝑒cos𝐸|||𝐸𝑓𝐸0,𝑎𝑒𝑒=𝑎5𝜇31𝑒2525𝑀4𝑒sin3𝐸+4sin21𝐸12𝑒sin3𝐸|||𝐸𝑓𝐸0,𝑎𝑒𝜔=𝑎𝜔𝑒=𝑎5𝜇31𝑒21/2𝑒54𝑒cos1𝐸+4𝑒23cos21𝐸+12𝑒cos3𝐸|||𝐸𝑓𝐸0,𝑎𝜔𝜔=𝑎5𝜇31𝑒2522𝑒2𝑀+𝑒54𝑒2sin1𝐸232𝑒2sin21𝐸+12𝑒sin3𝐸|||𝐸𝑓𝐸0,(7.3), 𝐸 and 𝑀=𝑛(𝑡𝜏) replacing 𝑛=𝜇/𝑎3, 𝜏, 𝑂0(𝑎0,𝑒0,𝜔0), and can be put in a suitable matrix form 𝑂𝑓(𝑎𝑓,𝑒𝑓,𝜔𝑓)where 𝑡𝑓 denotes the imposed changes on orbital elements (state variable), 𝑝0, 𝑡𝑓𝑡0=5, 𝑎 is the 𝑒 matrix of initial values of the adjoint variable and 𝜔 is 𝐹1 symmetric matrix. The adjoint variables are constant, and the matrix 𝑆1 is given by 𝑂 where 𝐻(𝑥,𝑦,𝜀) where 𝜀 is calculated by solving Kepler’s equation for 𝐻(𝑥,𝑦,𝜀)=𝐻0(𝑥,𝑦)+𝑘=1𝜀𝑘𝐻𝑘(𝑥,𝑦),(A.1), with 𝐻0(𝑥,𝑦). 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 [1517]; 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.

tab3
Table 3: Set of terminal orbits.
fig11
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 𝐹𝑥,𝑦=𝐹0𝑥,𝑦+𝑘=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 𝜀, 𝐻0𝑥,𝑦=𝐹0𝑥,𝑦,(A.5)This transformation is such that the new dynamical system has some advantages for the solution. Here, prime denotes the new variables, and 𝐻0,𝑆1+𝐻1=𝐹1,(A.6) and 𝐻0,𝑆2+12𝐻1+𝐹1,𝑆1+𝐻2=𝐹2,.(A.7) are also (𝑥,𝑦)-vectors. The new Hamiltonian 𝑚, resulting from this canonical transformation, is also developed as a power series of the small parameter 𝐻0,𝑆𝑚+Θ𝑚=𝐹𝑚,(A.8), Θ𝑚

The transformation equations are given explicitly by 𝐻0 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: 𝑐1,,𝑐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.

References

  1. J.-P. Marec, Optimal Space Trajectories, Elsevier, New York, NY, USA, 1979. View at Zentralblatt MATH
  2. G. Hori, “Theory of general perturbations with unspecified canonical variables,” Publications of the Astronomical Society of Japan, vol. 18, no. 4, pp. 287–296, 1966. View at Google Scholar
  3. C. Lánczos, The Variational Principles of Mechanics, Mathematical Expositions, no. 4, University of Toronto Press, Toronto, Canada, 4th edition, 1970. View at Zentralblatt MATH · View at MathSciNet
  4. T. N. Edelbaum, “Optimum power-limited orbit transfer in strong gravity fields,” AIAA Journal, vol. 3, pp. 921–925, 1965. View at Google Scholar · View at MathSciNet
  5. J.-P. Marec and N. X. Vinh, “Optimal low-thrust, limited power transfers between arbitrary elliptical orbits,” Acta Astronautica, vol. 4, no. 5-6, pp. 511–540, 1977. View at Publisher · View at Google Scholar
  6. L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, The Mathematical Theory of Optimal Processes, John Wiley & Sons, New York, NY, USA, 1962. View at Zentralblatt MATH · View at MathSciNet
  7. S. Da Silva Fernandes, “Generalized canonical systems—III. Space dynamics applications: solution of the coast-arc problem,” Acta Astronautica, vol. 32, no. 5, pp. 347–354, 1994. View at Publisher · View at Google Scholar
  8. R. H. Battin, An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Education Series, American Institute of Aeronautics and Astronautics, Washington, DC, USA, 1987. View at Zentralblatt MATH · View at MathSciNet
  9. T. N. Edelbaum, “An asymptotic solution for optimum power limited orbit transfer,” AIAA Journal, vol. 4, no. 8, pp. 1491–1494, 1966. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. G. E. Forsythe, M. A. Malcolm, and C. B. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, Englewood Cliffs, NJ, USA, 1977. View at Zentralblatt MATH · View at MathSciNet
  11. J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, vol. 12 of Texts in Applied Mathematics, Springer, New York, NY, USA, 3rd edition, 2002. View at Zentralblatt MATH · View at MathSciNet
  12. G. A. Bliss, Lectures on the Calculus of Variations, University of Chicago Press, Chicago, Ill, USA, 1946. View at Zentralblatt MATH · View at MathSciNet
  13. J.-P. Marec and N. X. Vinh, Étude Générale des Transferts Optimaux à Poussée Faible et Puissance Limitée entre Orbites Elliptiques Quelconques, ONERA, Châtillon, France, 1980.
  14. C. M. Haissig, K. D. Mease, and N. X. Vinh, “Minimum-fuel power-limited transfers between coplanar elliptical orbits,” Acta Astronautica, vol. 29, no. 1, pp. 1–15, 1993. View at Publisher · View at Google Scholar
  15. T. N. Edelbaum, “Optimum low-thrust rendezvous and station keeping,” AIAA Journal, vol. 2, no. 7, pp. 1196–1201, 1964. View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  16. J.-P. Marec, Transferts Optimaux Entre Orbites Elliptiques Proches, ONERA, Châtillon, France, 1967.
  17. S. Da Silva Fernandes, “Optimum low-thrust limited power transfers between neighbouring elliptic non-equatorial orbits in a non-central gravity field,” Acta Astronautica, vol. 35, no. 12, pp. 763–770, 1995. View at Publisher · View at Google Scholar
  18. S. Da Silva Fernandes, “Notes on Hori method for canonical systems,” Celestial Mechanics & Dynamical Astronomy, vol. 85, no. 1, pp. 67–77, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet