Analytical Solution for Optimal Low-Thrust Limited-Power Transfers Between Non-Coplanar Coaxial Orbits

ABSTRACT:

In this paper, an analytical solution for time-fixed optimal low-thrust limited-power transfers (no rendezvous) between elliptic coaxial non-coplanar orbits in an inverse-square force field is presented. Two particular classes of maneuvers are related to such transfers: maneuvers with change in the inclination of the orbital plane and maneuvers with change in the longitude of the ascending node. The optimization problem is formulated as a Mayer problem of optimal control with the state defined by semi-major axis, eccentricity, inclination or longitude of the ascending node, according to the class of maneuver considered, and a variable measuring the fuel consumption. After applying Pontryagin's maximum principle and determining the maximum Hamiltonian, short periodic terms are eliminated through an infinitesimal canonical transformation. The new maximum Hamiltonian resulting from this canonical transformation describes the extremal trajectories for long duration transfers. Closed-form analytical solution is then obtained through Hamilton-Jacobi theory. For long duration maneuvers, the existence of conjugate points is investigated through the Jacobi condition. Simplified solution is determined for transfers between close orbits. The analytical solution is compared to the numerical solution obtained by integration of the canonical system of differential equations describing the extremal trajectories for some sets of initial conditions. Results show a great agreement between these solutions for the class of maneuvers considered in the analysis. The solution of the two-point boundary value problem of going from an initial orbit to a final orbit, based on the analytical solution, is also discussed.

Keywords:
Low-thrust limited-power trajectories; Transfers between non-coplanar orbits; Optimal space trajectories

INTRODUCTION

In the last thirty years, important space missions have made use of low-thrust propulsion systems. The two pioneer missions were NASA-JPL Deep Space 1 and ESA's SMART-1. Deep Space 1 was the first interplanetary spacecraft to utilize Solar Electric Propulsion. It was developed by NASA in the New Millennium Program to testing new technologies for future Space and Earth Science Programs. It was launched on October 24, 1998, and its mission terminated on December 18, 2001, when its fuel supply exhausted. SMART-1 was the first of a series of ESA's Small Missions for Advanced Research in Technology. It was used to test solar electric propulsion and other deep-space technologies. It was launched on September 27, 2003, and its mission ended on September 3, 2006, when the spacecraft, in a planned maneuver, impacted the lunar surface. Interesting details about these space missions can be found in Rayman et al. (2000), Racca et al. (2002) and Camino et al. (2005).

Low-thrust electric propulsion systems are characterized by high specific impulse and low-thrust capability (the ratio between the maximum thrust acceleration and the gravity acceleration on the ground is small, between 10^-4 and 10^-2) and have their greatest benefits for high-energy planetary missions (Marec 1979; Racca 2003). Several researchers have obtained numerical and analytical solutions for several maneuvers involving specific initial and final orbits and specific thrust profiles (Gobetz 1964, 1965; Edelbaum 1964, 1965, 1966; Marec and Vinh 1977, 1980; Haissig et al. 1993; Geffroy and Epenoy 1997; Sukhanov 2000; Kiforenko 2005; Bonnard et al. 2006; da Silva Fernandes and das Chagas Carvalho 2008; Jamison and Coverstone, 2010; Jiang et al. 2012; Huang 2012; Quarta and Mengali 2013). In the most of the analytical studies, averaging techniques are applied and solutions of the averaged equations are obtained such that only secular behavior of the optimal solutions is discussed. Few works discuss the inclusion of periodic terms, which are, in general, considered only for transfers between close orbits (Gobetz 1964, 1965; Edelbaum 1965, 1966).

In two previous works (da Silva Fernandes and das Chagas Carvalho 2008; da Silva Fernandes et al. 2015), the authors presented complete analytical solutions based on canonical transformation, which include the short periodic terms, for the problem of optimal time-fixed low-thrust limited-power transfers between arbitrary elliptic coplanar orbits and between coplanar orbits with small eccentricities.

The main goal of this work is to extend the previous results and develop an approximated analytical solution, which includes short periodic terms, for the problem of optimal time-fixed low-thrust limited-power transfers (no rendezvous) between non-coplanar coaxial orbits in a Newtonian central gravity field. Two particular classes of maneuvers are related to such transfers: maneuvers with change in the inclination of the orbital plane and maneuvers with change in the longitude of the ascending node. The optimization problem is formulated as a Mayer problem of optimal control theory with semi-major axis, eccentricity, inclination or longitude of the ascending node, and a variable measuring the fuel consumption as state variables. Both maneuvers are described by a simplified set of the well-known Gauss equations of Celestial Mechanics. Pontryagin's maximum principle (Pontryagin et al. 1962) is applied to determine candidates for optimal trajectories, called extremals. After applying the set of necessary conditions provided by Pontryagin's maximum principle and determining the maximum Hamiltonian which describes the extremal trajectories, the short periodic terms are eliminated through an infinitesimal canonical transformation built using Hori method - a perturbation canonical method based on Lie series (Hori 1966). The new maximum Hamiltonian resulting from the infinitesimal canonical transformation describes the extremal trajectories for long duration transfers. A set of first integrals is derived for the canonical system of differential equations governed by the new Hamiltonian function. Closed-form analytical solutions are then obtained through Hamilton-Jacobi theory. The separation of variables technique is applied to solve the Hamilton-Jacobi equation. For long duration maneuvers, the existence of conjugate points is investigated through the Jacobi condition. 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 orbits a simplified solution is straightforwardly derived by linearizing the new Hamiltonian and the generating function obtained through the Hori method. This simplified solution is given by a linear system of algebraic equations in the initial value of the adjoint variables, such that the two-point boundary value problem can be solved by simple techniques. The analytical solutions are compared to the numerical solutions obtained by integration of the canonical system of differential equations describing the extremal trajectories for some sets of initial conditions.

It is noteworthy to mention that the results presented in this paper and obtained through canonical transformation theory are a complement of the works by Edelbaum (1965) and Marec and Vinh (1977, 1980).

FORMULATION OF THE OPTIMIZATION PROBLEM

A low-thrust limited-power propulsion system - LP system - is characterized by low-thrust acceleration level and high specific impulse. For such a system, the fuel consumption is described by the variable J defined as (Marec 1979):

(1) $\mathit{J}=\frac{1}{2}{\int }_{{\mathit{t}}_0}^{\mathit{t}}{\mathit{\gamma }}^2\mathit{dt},$

where: γ is the magnitude of the thrust acceleration vector γ, used as a control variable. The consumption variable J is a monotonic decreasing function of the mass m of the space vehicle,

$\mathit{J}={\mathit{P}}_{\max }\left(\frac{1}{\mathit{m}}-\frac{1}{{\mathit{m}}_0}\right)$

where: P_max is the maximum power and m₀ is the initial mass. The minimization of the final value of the fuel consumption, J_f, is equivalent to the maximization of m_f.

The motion of a space vehicle P, powered by a limited-power engine in an inverse-square force field, can be described by the well-known Gauss equations of Celestial Mechanics. These equations are given by Battin (1987):

(2) $\frac{\mathit{da}}{\mathit{dt}}=\frac{2}{\mathit{n}\sqrt{1-{\mathit{e}}^2}}\left[\mathit{e}\sin \mathit{fR}+\left(1+\mathit{e}\cos \mathit{f}\right)\mathit{S}\right],$

(3) $\frac{\mathit{de}}{\mathit{dt}}=\frac{\sqrt{1-{\mathit{e}}^2}}{\mathit{na}}\left[\sin \mathit{fR}+\left(\cos \mathit{E}+\cos \mathit{f}\right)\mathit{S}\right],$

(4) $\frac{\mathit{dI}}{\mathit{dt}}=\frac{1}{\mathit{na}\sqrt{1-{\mathit{e}}^2}}\left(\frac{\mathit{r}}{\mathit{a}}\right)\cos \left(\mathit{\omega }+\mathit{f}\right)\mathit{W},$

(5) $\frac{\mathit{d}\Omega }{\mathit{dt}}=\frac{1}{\mathit{na}\sqrt{1-{\mathit{e}}^2}\sin \mathit{I}}\left(\frac{\mathit{r}}{\mathit{a}}\right)\sin \left(\mathit{\omega }+\mathit{f}\right)\mathit{W},$

(6) $\frac{\mathit{d}\mathit{\omega }}{\mathit{dt}}=\frac{\sqrt{1-{\mathit{e}}^2}}{\mathit{nae}}\left[-\cos \mathit{f}\mathit{R}+\sin \mathit{f}\left(1+\frac{1}{1+\mathit{e}\cos \mathit{f}}\right)\mathit{S}\right]-\frac{\cos \mathit{I}}{\mathit{na}\sqrt{1-{\mathit{e}}^2}\sin \mathit{I}}\left(\frac{\mathit{r}}{\mathit{a}}\right)\sin \left(\mathit{\omega }+\mathit{f}\right)\mathit{W},$

(7) $\frac{\mathit{dM}}{\mathit{dt}}=\mathit{n}+\frac{1-{\mathit{e}}^2}{\mathit{nae}}\left[\left(\cos \mathit{f}-\frac{2\mathit{e}}{1+\mathit{e}\cos \mathit{f}}\right)\mathit{R}-\sin \mathit{f}\left(1+\frac{1}{1+\mathit{e}\cos \mathit{f}}\right)\mathit{S}\right],$

where: a, e, I, Ω, ω, M are the classical orbital elements: a is the semi-major axis, e is the eccentricity, I is the inclination of the orbital plane, Ω is the longitude of the ascending node, ω is the argument of periapsis, and M is the mean anomaly; and r is the radial distance. R, S and W are, respectively, radial, circumferential and normal component of the thrust acceleration. E is the eccentric anomaly, f is the true anomaly and n is the mean motion, n² a³ = µ, with µ denoting the gravitational parameter.

According to Marec and Vinh (1977, 1980), the optimal transfers between coaxial orbits are coaxial transfers, that is, the argument of periapsis does not change during the transfer, and two classes of maneuvers are related to the transfers between non-coplanar coaxial orbits: change in the inclination of the orbital plane and change in the longitude of the ascending node. For the first maneuver, ω = 0º or ω = 180º and Ω = 0º, and, for the second one, ω = 90º or 270º and I = 90º. A sketch of the geometric configuration of the initial orbit and the final orbit for these maneuvers is represented in Fig.1 (only "visible" portion of the orbits above the reference plane is represented).

Figure 1
Geometric configuration of the initial orbit and the final orbit: (a) represents the maneuver with change in the inclination of the orbital plane and (b) represents the maneuver with change in the longitude of ascending node.

Therefore, at time t, the state of a space vehicle P is defined by five variables: a, e, I or W, M and J. The inclination of the orbital plane or the longitude of the ascending node is used according to the class of maneuver considered in the analysis.

The optimization problem can be formulated as a Mayer problem of optimal control as follows:

It is proposed to transfer the space vehicle from the initial state (a₀, e₀, I₀ or W₀, M₀, 0) at time t₀ to the final state (a_f, e_f, I_f or W_f, M_f, J_f) at time t_f, such that the final consumption variable J_f is a minimum. The duration of the transfer t_f - t₀ is specified. For simple transfers (no rendezvous) the mean anomaly at time t_f is free.

So, taking into account the characterization of the two classes of maneuvers related to the transfers between non-coplanar coaxial orbits, described in the preceding paragraphs, one finds that the state equations are given by Gauss' equations for semi-major axis, eccentricity and mean anomaly, as defined by Eqs. 2, 3 and 7, respectively, and, by the following differential equations:

(8) $\frac{\mathit{d}\mathit{\theta }}{\mathit{dt}}=\pm \frac{1}{\mathit{na}\sqrt{1-{\mathit{e}}^2}}\left(\frac{\mathit{r}}{\mathit{a}}\right)\cos \mathit{f}\mathit{W},$

(9) $\frac{\mathit{dJ}}{\mathit{dt}}=\frac{1}{2}\left({\mathit{R}}^2+{\mathit{S}}^2+{\mathit{W}}^2\right).$

The variable θ is introduced to denote the inclination of the orbital plane or the longitude of the ascending node according to the maneuver considered. The upper (lower) sign refers to the maneuver with change in the inclination of the orbital plane with ω = 0º (ω = 180º) or maneuver with change in the longitude of the ascending node with ω = 90º (270º ). Note that Eq. 8 is derived straightforwardly form Eqs. 4 or 5 following the conditions described for each maneuver - maneuver with change in the inclination of the orbital plane or maneuver with change in the longitude of the ascending node - as described in a preceding paragraph.

According to the Pontryagin's maximum principle (Pontryagin et al. 1962), the adjoint variables (Lagrange multipliers associated to the constraints represented by state equations) p_a, p_e, p_θ, p_M and p_J are introduced and the Hamiltonian function H (a, e, θ, M, J, p_a, p_e, p_θ, p_M, p_J, R, S, W) is formed using the right-hand side of Eqs. 8 and 9,

(10) $\begin{array}{c}\mathit{H}={\mathit{np}}_{\mathit{M}}+\frac{1}{\mathit{na}\sqrt{1-{\mathit{e}}^2}}\left\{2\mathit{a}\left[\mathit{e}\sin \mathit{f}\mathit{R}+\left(1+\mathit{e}\cos \mathit{f}\right)\mathit{S}\right]{\mathit{p}}_{\mathit{a}}+\left(1-{\mathit{e}}^2\right)\left[\sin \mathit{f}\mathit{R}+\left(\cos \mathit{E}+\cos \mathit{f}\right)\mathit{S}\right]{\mathit{p}}_{\mathit{e}}\right.\hfill \\ +\frac{{\left(1-{\mathit{e}}^2\right)}^{3/2}}{\mathit{e}}\left[\left(\cos \mathit{f}-\frac{2\mathit{e}}{1+\mathit{e}\cos \mathit{f}}\right)\mathit{R}-\sin \mathit{f}\left(1+\frac{1}{1+\mathit{e}\cos \mathit{f}}\right)\mathit{S}\right]{\mathit{p}}_{\mathit{M}}+\left[\pm \left(\frac{\mathit{r}}{\mathit{a}}\right)\cos \mathit{f}\mathit{W}\right]\left.{\mathit{p}}_{\mathit{\theta }}\right\}+\frac{1}{2}\left({\mathit{R}}^2+{\mathit{S}}^2+{\mathit{W}}^2\right){\mathit{p}}_{\mathit{j}}.\hfill \end{array}$

The control variables R, S and W must be selected from the admissible controls such that the Hamiltonian function reaches its maximum along the optimal trajectory. Thus, the components of the optimal thrust acceleration are given by:

(11) ${\mathit{R}}^{*}=-\frac{1}{{\mathit{p}}_{\mathit{J}}\mathit{na}\sqrt{1-{\mathit{e}}^2}}\left[\left(2\mathit{ae}\sin \mathit{f}\right){\mathit{p}}_{\mathit{a}}+\left(1-{\mathit{e}}^2\right)\sin \mathit{f}{\mathit{p}}_{\mathit{e}}+\frac{{\left(1-{\mathit{e}}^2\right)}^{3/2}}{\mathit{e}}\left(\cos \mathit{f}-\frac{2\mathit{e}}{1+\mathit{e}\cos \mathit{f}}\right){\mathit{p}}_{\mathit{M}}\right],$

(12) ${\mathit{S}}^{*}=-\frac{1}{{\mathit{p}}_{\mathit{j}}\mathit{na}\sqrt{1-{\mathit{e}}^2}}\left[2\mathit{a}\left(1+\mathit{e}\cos \mathit{f}\right){\mathit{p}}_{\mathit{a}}+\left(1-{\mathit{e}}^2\right)\left(\cos \mathit{E}+\cos \mathit{f}\right){\mathit{p}}_{\mathit{e}}-\frac{{\left(1-{\mathit{e}}^2\right)}^{3/2}}{\mathit{e}}\sin \mathit{f}\left(1+\frac{1}{1+\mathit{e}\cos \mathit{f}}\right){\mathit{p}}_{\mathit{M}}\right],$

(13) ${\mathit{W}}^{*}=\mp \frac{1}{{\mathit{p}}_{\mathit{j}}\mathit{na}\sqrt{1-{\mathit{e}}^2}}\left(\frac{\mathit{r}}{\mathit{a}}\right)\cos \mathit{f}{\mathit{p}}_{\mathit{\theta }}.$

These equations can be simplified if it will be taken into account that the adjoint variable p_J is a first integral of the canonical system governed by the maximum Hamiltonian; that is, p_J is a constant whose value is obtained from the transversality condition, p_J (t_J) = -1. Accordingly, p_J (t_J) = -1. So, introducing this result into Eqs. 11 and 12, one finds the components of the optimal thrust acceleration.

Introducing Eqs. 11-13 into the Eq. 10, one finds:

(14) ${\mathit{H}}^{*}={\mathit{H}}_0+{\mathit{H}}_{\mathit{\gamma }}^{*},$

where H₀ denotes the undisturbed Hamiltonian and H*_γ is the part related to the optimal thrust acceleration. These functions are expressed by:

(15) $\begin{array}{c}{\mathit{H}}_0={\mathit{np}}_{\mathit{M}},\\ {\mathit{H}}_{\mathit{\gamma }}^{*}=\frac{1}{2{\mathit{n}}^2{\mathit{a}}^2\left(1-{\mathit{e}}^2\right)}\left\{\frac{1}{2}\left(1-\cos 2\mathit{f}\right){\left[2{\mathit{aep}}_{\mathit{a}}+\left(1-{\mathit{e}}^2\right){\mathit{p}}_{\mathit{e}}\right]}^2+4{\left(1-{\mathit{e}}^2\right)}^{3/2}\sin \mathit{f}\left(\frac{-2\mathit{e}}{1+\mathit{e}\cos \mathit{f}}+\cos \mathit{f}\right)\right.\end{array}$

(16) $\begin{array}{c}\left[{\mathit{ap}}_{\mathit{a}}{\mathit{p}}_{\mathit{M}}+\frac{\left(1-{\mathit{e}}^2\right)}{2\mathit{e}}{\mathit{p}}_{\mathit{e}}{\mathit{p}}_{\mathit{M}}\right]+\frac{{\left(1-{\mathit{e}}^2\right)}^3}{{\mathit{e}}^2}{\left(\frac{-2\mathit{e}}{1+\mathit{e}\cos \mathit{f}}+\cos \mathit{f}\right)}^2\mathit{p}_{\mathit{M}}{}^2+4{\mathit{a}}^2{\left(1-{\mathit{e}}^2\right)}^2{\left(\frac{\mathit{a}}{\mathit{r}}\right)}^2\mathit{p}_{\mathit{a}}{}^2+4\mathit{a}{\left(1-{\mathit{e}}^2\right)}^2\left(\frac{\mathit{a}}{\mathit{r}}\right)\\ \left(\cos \mathit{E}+\cos \mathit{f}\right){\mathit{p}}_{\mathit{a}}{\mathit{p}}_{\mathit{e}}+{\left(1-{\mathit{e}}^2\right)}^2{\left(\cos \mathit{E}+\cos \mathit{f}\right)}^2{\mathit{p}}_{\mathit{e}}^2-\frac{4\mathit{a}{\left(1-{\mathit{e}}^2\right)}^{5/2}}{\mathit{e}}\left(\frac{\mathit{a}}{\mathit{r}}\right)\sin \mathit{f}\left(1+\frac{1}{1+\mathit{e}\cos \mathit{f}}\right){\mathit{p}}_{\mathit{a}}{\mathit{p}}_{\mathit{M}}-\frac{2{\left(1-{\mathit{e}}^2\right)}^{5/2}}{\mathit{e}}\\ \left.\left(\cos \mathit{E}+\cos \mathit{f}\right)\sin \mathit{f}\left(1+\frac{1}{1+\mathit{e}\cos \mathit{f}}\right){\mathit{p}}_{\mathit{e}}{\mathit{p}}_{\mathit{M}}+{\left[-\frac{{\left(1-{\mathit{e}}^2\right)}^{3/2}}{\mathit{e}}\sin \mathit{f}\left(1+\frac{1}{1+\mathit{e}\cos \mathit{f}}\right){\mathit{p}}_{\mathit{M}}\right]}^2+{\left(\frac{\mathit{r}}{\mathit{a}}\right)}^2{\cos }^2\mathit{f}{\mathit{p}}_{\mathit{\theta }}^2\right\}.\end{array}$

It also be noted that the part of the maximum Hamiltonian related to the optimal thrust acceleration can be expressed in a compact form:

${\mathit{H}}_{\mathit{\gamma }}^{*}=\frac{1}{2}\mathit{\gamma }{*}^2.$

This result is used to compute the fuel consumption by simple quadrature of Eq. 9.

The problem of determining a first-order analytical solution of the system of differential equations governed by the maximum Hamiltonian H* is solved by means of the theory of canonical transformations as it will be described later.

FIRST-ORDER ANALYTICAL SOLUTION

A first-order analytical solution of the canonical system governed by the maximum Hamiltonian H* is derived by means of canonical transformation theory. Firstly, the short periodic terms are eliminated through an infinitesimal canonical transformation built using Hori method (Hori 1966). Hori method is a perturbation method based on Lie commutation theorem (Lie 1888) and it provides an explicit canonical transformation between old and new canonical variables, differently from the classical perturbation methods based on Hamilton-Jacobi theory, in which the transformation of variables involves new and old variables in a mixed form and requires the solution of an inversion problem to get the final form of the transformation of variables. We note that the new maximum Hamiltonian resulting from the infinitesimal canonical transformation describes the extremal trajectories for long duration transfers. A set of first integrals is then derived for the canonical system of differential equations governed by the new Hamiltonian function and a closed-form analytical solution is obtained through Hamilton-Jacobi theory. The separation of variables technique is applied to solve the Hamilton-Jacobi equation (Lanczos 1970).

ELIMINATION OF SHORT PERIODIC TERMS

In order to eliminate the short periodic terms from the maximum Hamiltonian, it is assumed that the undisturbed Hamiltonian H₀ is of zero order and the disturbing part H^*_γ is of the first-order in a small parameter defined by the magnitude of the thrust acceleration.

Consider an infinitesimal canonical transformation,

$\left(\mathit{a},\mathit{e},\mathit{\theta },\mathit{M},{\mathit{p}}_{\mathit{a}},{\mathit{p}}_{\mathit{e}},{\mathit{p}}_{\mathit{\theta }},{\mathit{p}}_{\mathit{M}}\right)\to \left({\mathit{a}}^{\prime },{\mathit{e}}^{\prime },{\mathit{\theta }}^{\prime },{\mathit{M}}^{\prime },{{\mathit{p}}^{\prime }}_{\mathit{a}},{{\mathit{p}}^{\prime }}_{\mathit{e}},{{\mathit{p}}^{\prime }}_{\mathit{\theta }},{{\mathit{p}}^{\prime }}_{\mathit{M}}\right).$

The new variables are designated by the prime.

According to the algorithm of Hori method, at order 0, one finds:

${\mathit{F}}_0={\mathit{n}}^{\prime }{{\mathit{p}}^{\prime }}_{\mathit{M}}.$

F₀ denotes the new undisturbed Hamiltonian.

In order to determine the higher order terms of the new Hamiltonian and the generating function of the canonical transformation, the algorithm proposed by da Silva Fernandes (2003) is applied. According to this algorithm, based on the method of variation of constants, the Poisson brackets involving the generating function and the new undisturbed Hamiltonian, that defines the general equation of the algorithm, is converted into a time partial differential equation. This procedure simplifies the determination of the new Hamiltonian and the generating function, since it eliminates the Hori auxiliary parameter.

So, consider the canonical system described by F₀:

$\begin{array}{cc}\frac{\mathit{d}{\mathit{a}}^{\prime }}{\mathit{dt}}=0& \frac{\mathit{d}{{\mathit{p}}^{\prime }}_{\mathit{a}}}{\mathit{dt}}=\frac{3}{2}\frac{{\mathit{n}}^{\prime }}{{\mathit{a}}^{\prime }}{{\mathit{p}}^{\prime }}_{\mathit{M}}\\ \frac{\mathit{d}{\mathit{e}}^{\prime }}{\mathit{dt}}=0& \frac{\mathit{d}{{\mathit{p}}^{\prime }}_{\mathit{e}}}{\mathit{dt}}=0\\ \frac{\mathit{d}{\mathit{\theta }}^{\prime }}{\mathit{dt}}=0& \frac{\mathit{d}{{\mathit{p}}^{\prime }}_{\mathit{\theta }}}{\mathit{dt}}=0,\\ \frac{\mathit{d}{\mathit{M}}^{\prime }}{\mathit{dt}}={\mathit{n}}^{\prime }& \frac{\mathit{d}{{\mathit{p}}^{\prime }}_{\mathit{M}}}{\mathit{dt}}0,\end{array}$

general solution of which is given by:

$\begin{array}{cc}{\mathit{a}}^{\prime }={{\mathit{a}}^{\prime }}_0& {{\mathit{p}}^{\prime }}_{\mathit{a}}={{\mathit{p}}^{\prime }}_{{\mathit{a}}_0}+\frac{3}{2}\frac{{\mathit{n}}^{\prime }\left(\mathit{t}-{\mathit{t}}_0\right)}{{\mathit{a}}^{\prime }}{{\mathit{p}}^{\prime }}_{\mathit{M}}\\ {\mathit{e}}^{\prime }={{\mathit{e}}^{\prime }}_0& {{\mathit{p}}^{\prime }}_{\mathit{e}}={{\mathit{p}}^{\prime }}_{{\mathit{e}}_0}\\ {\mathit{\theta }}^{\prime }={{\mathit{\theta }}^{\prime }}_0& {{\mathit{p}}^{\prime }}_{\mathit{\theta }}={{\mathit{p}}^{\prime }}_{{\mathit{\theta }}_0}\\ {\mathit{M}}^{\prime }={{\mathit{M}}^{\prime }}_0+{\mathit{n}}^{\prime }\left(\mathit{t}-{\mathit{t}}_0\right)& {{\mathit{p}}^{\prime }}_{\mathit{M}}={{\mathit{p}}^{\prime }}_{{\mathit{M}}_0}.\end{array}$

The subscript 0 denotes the constants of integration.

Introducing this general solution into the equation of order 1 of the algorithm of Hori method, one finds:

$\frac{\partial {\mathit{S}}_1}{\partial \mathit{t}}={\mathit{H}}_{\mathit{\gamma }}^{*}-{\mathit{F}}_1,$

with the functions H* _γ written in terms of the constants of integration.

Following the proposed algorithm, one finds:

(17) ${\mathit{F}}_1=\frac{{\mathit{a}}^{\prime }}{2\mathit{\mu }}\left\{4{{\mathit{a}}^{\prime }}^2{{\mathit{p}}^{\prime }}_{\mathit{a}}^2+\frac{5}{2}\left(1-{{\mathit{e}}^{\prime }}^2\right){{\mathit{p}}^{\prime }}_{\mathit{e}}^2+\frac{\left(1+4{{\mathit{e}}^{\prime }}^2\right)}{2\left(1-{{\mathit{e}}^{\prime }}^2\right)}{{\mathit{p}}^{\prime }}_{\mathit{\theta }}^2\right\},$

(18) $\begin{array}{c}{\mathit{S}}_1=\frac{1}{2}\sqrt{\frac{{{\mathit{a}}^{\prime }}^5}{{\mathit{\mu }}^3}}\left\{8{\mathit{e}}^{\prime }\right.\sin {\mathit{E}}^{\prime }{{\mathit{a}}^{\prime }}^2{{\mathit{p}}^{\prime }}_{\mathit{a}}^2+8\left(1-{{\mathit{e}}^{\prime }}^2\right)\sin {\mathit{E}}^{\prime }{\mathit{a}}^{\prime }{{\mathit{p}}^{\prime }}_{\mathit{a}}{{\mathit{p}}^{\prime }}_{\mathit{e}}+\left(1-{{\mathit{e}}^{\prime }}^2\right)\left[-\frac{5}{4}{\mathit{e}}^{\prime }\sin {\mathit{E}}^{\prime }+\frac{3}{4}\sin 2{\mathit{E}}^{\prime }-\frac{1}{12}{\mathit{e}}^{\prime }\sin 3{\mathit{E}}^{\prime }\right]{{\mathit{p}}^{\prime }}_{\mathit{e}}^2\hfill \\ +{\left(1-{{\mathit{e}}^{\prime }}^2\right)}^{-1}\left.\left[\left(-\frac{9}{4}{\mathit{e}}^{\prime }+{{\mathit{e}}^{\prime }}^3\right)\sin {\mathit{E}}^{\prime }+\left(\frac{1}{4}+\frac{1}{2}{{\mathit{e}}^{\prime }}^2\right)\sin 2{\mathit{E}}^{\prime }-\frac{1}{12}{\mathit{e}}^{\prime }\sin 3{\mathit{E}}^{\prime }\right]{{\mathit{p}}^{\prime }}_{\mathit{\theta }}^2\right\}.\hfill \end{array}$

Terms factored by p'_M have been omitted in equations above, since only transfers (no rendezvous) are considered. Note that p'_M is a first integral of the average canonical system; that is, p'_M is constant and its value is computed from the transversality condition which gives p'_M (t_f)= 0 for simple transfers (the final position of the space vehicle is free). This result simplifies the expressions of the new Hamiltonian and the generating function as given by equations above.

A SET OF FIRST INTEGRALS OF THE AVERAGED CANONICAL SYSTEM

Consider the Mathieu transformation (Lanczos 1970) defined by the following equations:

(19) $\begin{array}{cc}{\mathit{a}}^{\prime }={\mathit{a}}^{″}& {{\mathit{p}}^{\prime }}_{\mathit{a}}={{\mathit{p}}^{″}}_{\mathit{a}}\end{array}$

(20) $\begin{array}{cc}{\mathit{e}}^{\prime }=\sin \mathit{\phi }& {{\mathit{p}}^{\prime }}_{\mathit{e}}=\frac{{{\mathit{p}}^{″}}_{\mathit{\phi }}}{\cos \mathit{\phi }}\end{array}$

(21) $\begin{array}{cc}{\mathit{\theta }}^{\prime }={\mathit{\theta }}^{″}& {{\mathit{p}}^{\prime }}_{\mathit{\theta }}={{\mathit{p}}^{″}}_{\mathit{\theta }}.\end{array}$

The Hamiltonian function F₁ is invariant with respect to this transformation,

${\mathit{F}}^{″}=\frac{{\mathit{a}}^{″}}{2\mathit{\mu }}\left\{4{{\mathit{a}}^{″}}^2{{\mathit{p}}^{″}}_{\mathit{a}}^2+\frac{5}{2}{{\mathit{p}}^{″}}_{\mathit{\phi }}^2+\frac{1}{2}\left(1+5{\tan }^2\mathit{\phi }\right){{\mathit{p}}^{″}}_{\mathit{\theta }}^2\right\}.$

The averaged canonical system governed by the Hamiltonian function F'' has a set of first integrals, which can be derived straightforwardly from the differential equations as described by Pines (1964) and Edelbaum and Pines (1970).

The first integrals are given by the following equations:

(22) $\frac{{\mathit{a}}^{″}}{2\mathit{\mu }}\left\{4{{\mathit{a}}^{″}}^2{{\mathit{p}}^{″}}_{\mathit{a}}^2+\frac{5}{2}{{\mathit{p}}^{″}}_{\mathit{\phi }}^2+\frac{1}{2}\left(1+5{\tan }^2\mathit{\phi }\right){{\mathit{p}}^{″}}_{\mathit{\theta }}^2\right\}=E,$

(23) ${{\mathit{p}}^{″}}_{\mathit{\phi }}^2+{{\mathit{p}}^{″}}_{\mathit{\theta }}^2{\tan }^2\mathit{\phi }={\mathit{C}}_1^2,$

(24) ${{\mathit{p}}^{″}}_{\mathit{\theta }}={\mathit{C}}_2,$

E, C_i, i = 1,2 are constants.

CLOSED-FORM SOLUTION OF THE CANONICAL SYSTEM GOVERNED BY HAMILTONIAN FUNCTION F''

Consider the canonical transformation,

$\left({\mathit{a}}^{″},\mathit{\phi },{\mathit{\theta }}^{″},{{\mathit{p}}^{″}}_{\mathit{a}},{{\mathit{p}}^{″}}_{\mathit{\phi }},{{\mathit{p}}^{″}}_{\mathit{\theta }}\right)\stackrel{\mathit{W}}{\to }\left({\mathit{C}}_1,{\mathit{C}}_2,E,{\mathit{p}}_{{\mathit{c}}_1},{\mathit{p}}_{{\mathit{c}}_2},{\mathit{p}}_E\right),$

defined by a generating function W such that the constants C₁, C₂ and E become the new generalized coordinates.

Since the new Hamiltonian function F'' is a quadratic form in the adjoint variables or conjugate momenta, the separation of variables technique will be applied for solving the Hamilton-Jacobi equation (Lanczos 1970).

Consider the transformation equations:

(25) $\begin{array}{cc}{{\mathit{p}}^{″}}_{\mathit{a}}=\frac{\partial \mathit{W}}{\partial {\mathit{a}}^{″}}& {\mathit{p}}_{{\mathit{C}}_1}\end{array}=-\frac{\partial \mathit{W}}{\partial {\mathit{C}}_1},$

(26) $\begin{array}{cc}{\mathit{p}}_{\varphi }=\frac{\partial \mathit{W}}{\partial \varphi }& {\mathit{p}}_{{\mathit{C}}_2}\end{array}=-\frac{\partial \mathit{W}}{\partial {\mathit{C}}_2},$

(27) $\begin{array}{cc}{{\mathit{p}}^{″}}_{\mathit{\theta }}=\frac{\partial \mathit{W}}{\partial {\mathit{\theta }}^{″}}& {\mathit{p}}_E=-\frac{\partial \mathit{W}}{\partial E},\end{array}$

where: $\mathit{W}=\mathit{W}\left({\mathit{a}}^{″},\mathit{\phi },{\mathit{\theta }}^{″},{\mathit{C}}_1,{\mathit{C}}_2,E\right)$ is the generating function.

The corresponding Hamilton-Jacobi equation is then given by:

(28) $\frac{{\mathit{a}}^{″}}{2\mathit{\mu }}\left\{4{{\mathit{a}}^{″}}^2{\left(\frac{\partial \mathit{W}}{\partial {\mathit{a}}^{″}}\right)}^2+\frac{5}{2}{\left(\frac{\partial \mathit{W}}{\partial \mathit{\phi }}\right)}^2+\frac{1}{2}\left(1+5{\tan }^2\mathit{\phi }\right)\left(\frac{\partial \mathit{W}}{\partial {\mathit{\theta }}^{″}}\right)\right\}=E.$

Following the separation of variables technique (Lanczos 1970), it is assumed that the generating function W is equal to a sum of functions, each of which depends on a single old variable, that is:

(29) $\mathit{W}\left({\mathit{a}}^{″},\mathit{\phi },{\mathit{\theta }}^{″},{\mathit{C}}_1,{\mathit{C}}_2,E\right)={\mathit{W}}_1\left({\mathit{a}}^{″},{\mathit{C}}_1,{\mathit{C}}_2,E\right)+{\mathit{W}}_2\left(\mathit{\phi },{\mathit{C}}_1,{\mathit{C}}_2,E\right)+{\mathit{W}}_3\left({\mathit{\theta }}^{″},{\mathit{C}}_1,{\mathit{C}}_2,E\right).$

Therefore, from Eqs. 22-29, it follows that

$\begin{array}{c}{\left(\frac{\partial {\mathit{W}}_2}{{\partial }_{\mathit{\phi }}}\right)}^2+{\left(\frac{\partial {\mathit{W}}_3}{\partial {\mathit{\theta }}^{″}}\right)}^2{\tan }^2\mathit{\phi }={\mathit{C}}_1^2,\\ \frac{\partial {\mathit{W}}_3}{\partial {\mathit{\theta }}^{″}}={\mathit{C}}_2,\\ \frac{{\mathit{a}}^{″}}{2\mathit{\mu }}\left\{4{{\mathit{a}}^{″}}^2{\left(\frac{\partial {\mathit{W}}_1}{\partial {\mathit{a}}^{″}}\right)}^2+\frac{5}{2}{\left(\frac{\partial {\mathit{W}}_2}{\partial \mathit{\phi }}\right)}^2+\frac{1}{2}\left(1+5{\tan }^2\mathit{\phi }\right){\left(\frac{\partial {\mathit{W}}_3}{\partial {\mathit{\theta }}^{″}}\right)}^2\right\}=E.\end{array}$

A complete solution of these equations is given by:

(30) ${\mathit{W}}_1=-\sqrt{\frac{5{\mathit{C}}^2}{2}}\left\{{\left[\frac{4\mathit{\mu }E}{5{\mathit{C}}^2{\mathit{a}}^{″}}-1\right]}^{1/2}-{\tan }^{-1}{\left[\frac{4\mathit{\mu }E}{5\mathit{C}{}^2{\mathit{a}}^{″}}-1\right]}^{1/2}\right\},$

(31) ${\mathit{W}}_2=\mathbf{\int }\sqrt{{\mathit{C}}_1^2-{\mathit{C}}_2^2{\tan }^2\mathit{\phi }}\mathit{d}\mathit{\phi },$

(32) ${\mathit{W}}_3={\mathit{C}}_2{\mathit{\theta }}^{″},$

with $5{\mathit{C}}_1^2+{\mathit{C}}_2^2=5{\mathit{C}}^2$. We note that W₂ is given as indefinite integral, since only its partial derivatives are needed (see Eqs. 25-27), as shown in the next paragraphs.

Now, consider the differential equations for the conjugate momenta of the canonical system governed by the new Hamiltonian F'' = E. These equations are:

$\begin{array}{ccc}\frac{{\mathit{dp}}_{{\mathit{C}}_1}}{\mathit{dt}}=0& \frac{{\mathit{dp}}_{{\mathit{C}}_2}}{\mathit{dt}}=0& \frac{{\mathit{dp}}_E}{\mathit{dt}}=-1,\end{array}$

whose solution is very simple:

(33) $\begin{array}{ccc}{\mathit{p}}_{{\mathit{C}}_1}={\mathit{\alpha }}_1& {\mathit{p}}_{{\mathit{c}}_2}={\mathit{\alpha }}_2& {\mathit{p}}_E={\mathit{\alpha }}_3-\mathit{t},\end{array}$

where: α_i, i = 1, 2, 3 are arbitrary constants of integration.

Introducing the generating function W, defined by Eqs. 29-32, into the transformation equations, defined by Eqs. 25-27, and taking into account the general solution defined by Eqs. 33 for the conjugate momenta, one finds:

(34) ${{\mathit{p}}^{″}}_{\mathit{a}}=\frac{1}{2\sqrt{2}{{\mathit{a}}^{″}}^2}{\left[4\mathit{\mu }E{\mathit{a}}^{″}-5{\mathit{C}}^2{{\mathit{a}}^{″}}^2\right]}^{1/2},$

(35) ${\mathit{p}}_{\mathit{\phi }}=\sqrt{{\mathit{C}}_1^2-{\mathit{C}}_2^2{\tan }^2\mathit{\phi },}$

(36) ${{\mathit{p}}^{″}}_{\mathit{\theta }}={\mathit{C}}_2,$

(37) ${\mathit{\alpha }}_1=-\frac{{\mathit{C}}_1}{\mathit{C}}\sqrt{\frac{5}{2}}{\tan }^{-1}\left({\left(\frac{4\mathit{\mu }E}{5{\mathit{C}}^2{\mathit{a}}^{″}}-1\right)}^{1/2}\right)-\frac{{\mathit{C}}_1}{\sqrt{{\mathit{C}}_1^2+{\mathit{C}}_2^2}}{\tan }^{-1}\left(\frac{\sqrt{{\mathit{C}}_1^2+{\mathit{C}}_2^2}\tan \mathit{\phi }}{\sqrt{{\mathit{C}}_1^2-{\mathit{C}}_2^2{\tan }^2\mathit{\phi }}}\right),$

(38) ${\mathit{\alpha }}_2=-\frac{1}{5}\sqrt{\frac{5}{2}}\frac{{\mathit{C}}_2}{\mathit{C}}{\tan }^{-1}\left({\left(\frac{4\mathit{\mu }E}{5{\mathit{C}}^2{\mathit{a}}^{″}}-1\right)}^{1/2}\right)+{\tan }^{-1}\left(\frac{{\mathit{C}}_2\tan \mathit{\phi }}{\sqrt{{\mathit{C}}_1^2-{\mathit{C}}_2^2{\tan }^2\mathit{\phi }}}\right)-\frac{{\mathit{C}}_2}{\sqrt{{\mathit{C}}_1^2+{\mathit{C}}_2^2}}{\tan }^{-1}\left(\frac{\sqrt{{\mathit{C}}_1^2+{\mathit{C}}_2^2}\tan \mathit{\phi }}{\sqrt{{\mathit{C}}_1^2-{\mathit{C}}_2^2{\tan }^2\mathit{\phi }}}\right)-\mathit{\theta },$

(39) $\mathit{t}-{\mathit{\alpha }}_3=-\frac{1}{2\sqrt{2}E{\mathit{a}}^{″}}{\left[4\mathit{\mu }E{\mathit{a}}^{″}-5{\mathit{C}}^2{{\mathit{a}}^{″}}^2\right]}^{1/2}.$

From the above equations one finds, after some simplifications, the solution of the canonical system governed by the Hamiltonian F'' for a given set of initial conditions:

(40) ${\mathit{a}}^{″}\left(\mathit{t}\right)=\frac{{{\mathit{a}}^{″}}_0}{1+\frac{4{{\mathit{a}}^{″}}_0}{\mathit{\mu }}\left(\frac{1}{2}E{\mathit{t}}^2-{{\mathit{a}}^{″}}_{\mathit{a}}{{\mathit{p}}^{″}}_{{\mathit{a}}_0}\mathit{t}\right)},$

(41) ${\mathit{a}}^{″}{\sin }^2{\mathit{k}}_0={{\mathit{a}}_0}^{″}{\sin }^2\left(\sqrt{2}\mathit{\psi }+{\mathit{k}}_0\right),$

(42) $\mathit{\psi }={\frac{\mathit{C}}{\sqrt{5\left({\mathit{C}}_1^2+{\mathit{C}}_2^2\right)}}}^{\left(\mathit{\tau }-{\mathit{\tau }}_0\right),}$

(43) $\sin \mathit{\phi }=\sin {\mathit{k}}_1\sin \mathit{\tau },$

(44) ${\mathit{\theta }}^{″}={\mathit{k}}_2+{\tan }^{-1}\left(\tan \mathit{\tau }/\sec {\mathit{k}}_1\right)-\frac{4}{5}\mathit{\tau }\cos {\mathit{k}}_1,$

(45) ${{\mathit{p}}^{″}}_{\mathit{a}}^2={\left(\frac{{{\mathit{a}}_0}^{″}}{{\mathit{a}}^{″}}\right)}^3{{\mathit{p}}^{″}}_{{\mathit{a}}_0}^2+\frac{1}{8}{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2\left(5{\sec }^2{\mathit{k}}_1-4\right)\left(\frac{{{\mathit{a}}_0}^{″}}{{{\mathit{a}}^{″}}^3}-\frac{1}{{{\mathit{a}}^{″}}^2}\right),$

(46) ${{\mathit{p}}^{″}}_{\mathit{\phi }}^2={{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2\left({\sec }^2{\mathit{k}}_1-{\sec }^2\mathit{\phi }\right),$

(47) ${{\mathit{p}}^{″}}_{\mathit{\theta }}={{\mathit{p}}^{″}}_{{\mathit{\theta }}_0},$

with the auxiliary constants k₀, k₁ and k₂ defined as functions of the initial value of the adjoint variables by:

(48) ${\csc }^2{\mathit{k}}_0=\frac{8{\left({{\mathit{a}}^{″}}_0{{\mathit{p}}^{″}}_{{\mathit{a}}_0}\right)}^2+{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2\left(5{\sec }^2{\mathit{k}}_1-4\right)}{{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2\left(5{\sec }^2{\mathit{k}}_1-4\right)},$

(49) ${\sec }^2{\mathit{k}}_1=\frac{{{\mathit{p}}^{″}}_{{\mathit{\phi }}_0}^2+{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2{\sec }^2{\mathit{\phi }}_0}{{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2},$

(50) ${\mathit{k}}_2={{\mathit{\theta }}^{″}}_0-+{\tan }^{-1}\left(\tan {\mathit{\tau }}_0/\sec {\mathit{k}}_1\right)+\frac{4}{5}{\mathit{\tau }}_0\cos {\mathit{k}}_1.$

Combining Eqs. 34 and 39, one finds:

${\mathit{a}}^{″}{{\mathit{p}}^{″}}_{\mathit{a}}={{\mathit{a}}^{″}}_0{{\mathit{p}}^{″}}_{{\mathit{a}}_0}-\mathit{Et}.$

This equation can also be used to evaluate the adjoint variable p''_a along the trajectory, as well as Eq. 45.

The constants C, C₁, C₂ and E can also be written as functions of the initial value of the adjoint variables:

(51) ${\mathit{C}}^2=\frac{1}{5}{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2\left(5{\sec }^2{\mathit{k}}_1-4\right),$

(52) ${\mathit{C}}_1^2={{\mathit{p}}^{″}}_{{\mathit{\phi }}_0}^2+{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2{\tan }^2{\mathit{\phi }}_0,$

(53) ${\mathit{C}}_2={{\mathit{p}}^{″}}_{{\mathit{\theta }}_0},$

(54) $4\mathit{\mu }E={{\mathit{a}}^{″}}_0\left(8{\left({{\mathit{a}}^{″}}_0{{\mathit{p}}^{″}}_{{\mathit{a}}_0}\right)}^2+{{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2\left(5{\sec }^2{\mathit{k}}_1-4\right)\right).$

In all equations above, the initial conditions are defined by a'' (0) = a''₀, e'' (0) = sinφ₀ and θ''(0) = θ''₀, and τ₀ is obtained from the Eq. 43, that is, sinφ₀ = sink₁ sinτ₀.

An alternative expression for variable ψ can be determined from Eqs. 42, 49, 51, 52 and 53, as it follows:

(55) $\mathit{\psi }=\frac{1}{5}\sqrt{1+4{\sin }^2{\mathit{k}}_1}\left(\mathit{\tau }-{\mathit{\tau }}_0\right).$

This equation will be useful in the analysis of conjugate point for long duration transfers, discussed later. Note that ψ does not depend on k₀.

FIRST-ORDER SOLUTION OF THE CANONICAL SYSTEM GOVERNED BY HAMILTONIAN FUNCTION H*

Following the algorithm of Hori method (Hori 1966), the partial derivatives of the generating function S₁ with respect to the adjoint variables must be determined in order to obtain a first-order solution of the canonical system governed by Hamiltonian functionH*. Computing these partial derivatives and applying the initial conditions, one finds:

(56) $\mathit{a}\left(\mathit{t}\right)={\mathit{a}}^{\prime }\left(\mathit{t}\right)+\sqrt{\frac{{{\mathit{a}}^{\prime }}^5}{{\mathit{\mu }}^3}}{\left[8{\mathit{e}}^{\prime }\sin {\mathit{E}}^{\prime }{{\mathit{a}}^{\prime }}^2{{\mathit{p}}^{\prime }}_{\mathit{a}}+4\left(1-{{\mathit{e}}^{\prime }}^2\right)\sin {\mathit{E}}^{\prime }{\mathit{a}}^{\prime }{{\mathit{p}}^{\prime }}_{\mathit{e}}\right]}_{{\mathit{t}}_0}^{\mathit{t}},$

(57) $\mathit{e}\left(\mathit{t}\right)={\mathit{e}}^{\prime }\left(\mathit{t}\right)+\sqrt{\frac{{{\mathit{a}}^{\prime }}^5}{{\mathit{\mu }}^3}}{\left[4\left(1-{{\mathit{e}}^{\prime }}^2\right)\sin {\mathit{E}}^{\prime }{\mathit{a}}^{\prime }{{\mathit{p}}^{\prime }}_{\mathit{a}}+\left(1-{{\mathit{e}}^{\prime }}^2\right)\left[-\frac{5}{4}{\mathit{e}}^{\prime }\sin {\mathit{E}}^{\prime }+\frac{3}{4}\sin 2{\mathit{E}}^{\prime }-\frac{1}{12}{\mathit{e}}^{\prime }\sin 3{\mathit{E}}^{\prime }\right]{{\mathit{p}}^{\prime }}_{\mathit{e}}\right]}_{{\mathit{t}}_0}^{\mathit{t}},$

(58) $\mathit{\theta }\left(\mathit{t}\right)={\mathit{\theta }}^{\prime }\left(\mathit{t}\right)+\sqrt{\frac{{{\mathit{a}}^{\prime }}^5}{{\mathit{\mu }}^3}}{\left(1-{{\mathit{e}}^{\prime }}^2\right)}^{-1}\left[\left(-\frac{9}{4}{\mathit{e}}^{\prime }+{{\mathit{e}}^{\prime }}^3\right)\sin {\mathit{E}}^{\prime }+\left(\frac{1}{4}+\frac{1}{2}{{\mathit{e}}^{\prime }}^2\right)\sin 2{\mathit{E}}^{\prime }-\frac{1}{12}{\mathit{e}}^{\prime }\sin 3{\mathit{E}}^{\prime }\right]{\left.{{\mathit{p}}^{\prime }}_{\mathit{\theta }}\right]}_{{\mathit{t}}_0}^{\mathit{t}},$

with a', e', ..., p'_θ previously defined through the Mathieu transformation given by Eqs. 19-21. The eccentric anomaly E' is computed from the well-known Kepler's equation with the mean anomaly given by:

(59) ${\mathit{M}}^{\prime }\left(\mathit{t}\right)={\mathit{M}}^{\prime }\left({\mathit{t}}_0\right)+{\int }_{{\mathit{t}}_0}^{\mathit{t}}{\mathit{ndt}}^{\prime }.$

An approximate expression for the fuel consumption along the extremal trajectory is obtained by simple quadrature of Eq. 9 and can be put in a compact form as (da Silva Fernandes and das Chagas Carvalho 2008):

(60) $\mathit{J}=E\left(\mathit{t}-{\mathit{t}}_0\right)+\Delta {\mathit{S}}_1,$

with $\Delta S_1=S_1\left(M^{\prime }\right)-S_1\left({{\mathit{M}}^{\prime }}_0\right)$, with S₁ given by Eq. 18. Note that Eq. 60 includes the periodic terms through the generating function S₁.

LONG DURATION TRANSFERS

In what follows is presented a discussion about the existence of conjugate points through the analysis of Jacobi condition, and a description of an algorithm for solving the two-point boundary value problem of going form an initial orbit to a final orbit for long duration transfers. Such transfers are described by the general solution of the canonical system governed by the Hamiltonian F''

EXISTENCE OF CONJUGATE POINT

Consider Eqs. 41-44, which define a two-parameter family of extremals in the phase space (a'', φ'', θ'') for a given initial phase point (a₀'', φ₀, θ₀'') corresponding to an initial orbit. By eliminating the auxiliary variables τ and ψ (see Eq. 55), α =a''/a₀'' and θ'' can be written as explicit functions of φ, that is, α = α(φ, φ₀, k₀, k₁) and θ'' (φ, φ₀, θ₀'', k₁). The conjugate points to the phase point (a₀'', φ₀, θ₀'') are determined through the roots of the equation (Bliss 1946):

$\frac{\partial \mathit{\alpha }}{\partial {\mathit{k}}_0}\frac{\partial {\mathit{\theta }}^{″}}{\partial {\mathit{k}}_1}-\frac{\partial \mathit{\alpha }}{\partial {\mathit{k}}_1}\frac{\partial {\mathit{\theta }}^{″}}{\partial {\mathit{k}}_0}=0.$

Since θ'' does not depend on k₀, ∂θ''/∂k₀ = 0.

On the other hand, from Eq. 41 and from the remark about Eq. 55, one finds:

$\frac{\partial \mathit{\alpha }}{\partial {\mathit{k}}_0}=-\frac{2\sqrt{\mathit{\alpha }}}{\sin {\mathit{k}}_0}{\left(\mathit{\alpha }-2\sqrt{\mathit{\alpha }}\cos \left(\sqrt{2}\mathit{\psi }\right)+1\right)}^{1/2}\cancel{=}0.$

Thus, the problem of determining the conjugate points reduces to the analysis of the roots of the following equation:

$\frac{\partial {\mathit{\theta }}^{″}}{\partial {\mathit{k}}_1}=0.$

From Eqs. 43 and 44, it follows that

(61) ${\mathit{\theta }}^{″}={{\mathit{\theta }}^{″}}_0-\frac{4}{5}\cos \mathit{k}{}_1\left\{{\sin }^{-1}\left(\frac{\sin \mathit{\phi }}{\sin {\mathit{k}}_1}\right)-{\sin }^{-1}\frac{\sin {\mathit{\phi }}_0}{\sin {\mathit{k}}_1}\right\}+{\sin }^{-1}\left(\frac{\tan \mathit{\phi }}{\tan {\mathit{k}}_1}\right)-{\sin }^{-1}\left(\frac{\tan {\mathit{\phi }}_0}{\tan {\mathit{k}}_1}\right).$

The partial derivative ∂θ'' = ∂k₁ is then given by

(62) $\frac{\partial {\mathit{\theta }}^{″}}{\partial {\mathit{k}}_1}=-\frac{4}{5}\left(\mathit{\tau }-{\mathit{\tau }}_0\right)\sin {\mathit{k}}_1+\frac{4}{5}{\cot }^2{\mathit{k}}_1\left[-\frac{\sin {\mathit{\phi }}_0}{\cos {\mathit{\tau }}_0}+\frac{\sin \mathit{\phi }}{\cos \mathit{\tau }}\right]+\frac{{\sec }^2{\mathit{k}}_1}{\tan {\mathit{k}}_1}\left[\frac{\tan {\mathit{\phi }}_0}{\sqrt{{\tan }^2{\mathit{k}}_1-{\tan }^2{\mathit{\phi }}_0}}-\frac{\tan \mathit{\phi }}{\sqrt{{\tan }^2{\mathit{k}}_1-{\tan }^2{\mathit{\phi }}_0}}\right].$

The auxiliary variable τ is reintroduced to simplify the expression.

Therefore, after some simplifications, one finds that the conjugate points are given by the roots of following equation:

(63) $\left(1+4{\sin }^2{\mathit{k}}_1\right)\sin \left(\mathit{\tau }-{\mathit{\tau }}_0\right)+4\left(\mathit{\tau }-{\mathit{\tau }}_0\right){\sin }^2{\mathit{k}}_1\cos \mathit{\tau }\cos {\mathit{\tau }}_0=0.$

If a conjugate point τ* exists, or, equivalently, φ* (see Eq. 43 connecting τ and φ), then the extremal does not yield a local minimum. This is the necessary Jacobi condition in the Calculus of Variations for an extremum (Bliss 1946).

SOLUTION OF THE TWO-POINT BOUNDARY VALUE PROBLEM

In what follows, the solution of the two-point boundary value problem of going from an initial orbit to a final orbit for long duration transfers is discussed.

To describe this solution, a new consumption variable is introduced (Edelbaum 1965). This new variable is defined by:

$\mathit{u}=\sqrt{2E\mathit{t}.}$

Accordingly, Eqs. 40 and 41 can be put in the form:

(64) $\frac{{\mathit{a}}^{″}}{{{\mathit{a}}^{″}}_0}={\left[1-2\left(\frac{\mathit{u}}{{\mathit{v}}_0}\right)\cos {\mathit{k}}_0+{\left(\frac{\mathit{u}}{{\mathit{v}}_0}\right)}^2\right]}^{-1},$

(65) $\left(\frac{\mathit{u}}{{\mathit{v}}_0}\right)=\frac{\sin \left(\sqrt{2\mathit{\psi }}\right)}{\sin \left(\sqrt{2\mathit{\psi }+{\mathit{k}}_0}\right)},$

where ${\mathit{v}}_0=\sqrt{\mathit{\mu }/{{\mathit{a}}^{″}}_0}$. These equations are the same ones obtained for transfers between coplanar orbits (da Silva Fernandes and das Chagas Carvalho 2008).

By eliminating k₀ from the above equations, one finds:

(66) $\left(\frac{\mathit{u}}{{\mathit{v}}_0}\right)=\sqrt{1-2\sqrt{\frac{{{\mathit{a}}^{″}}_0}{{\mathit{a}}^{″}}}\cos \left(\sqrt{2\mathit{\psi }}\right)+\left(\frac{{{\mathit{a}}^{″}}_0}{{\mathit{a}}^{″}}\right)}.$

On the other hand, J = Et, thus:

(67) $\mathit{J}=\frac{{\mathit{v}}_0^2}{2\mathit{t}}\left[1-2\sqrt{\frac{{{\mathit{a}}^{″}}_0}{{\mathit{a}}^{″}}}\cos \left(\sqrt{2}\mathit{\psi }\right)+\left(\frac{{{\mathit{a}}^{″}}_0}{{\mathit{a}}^{″}}\right)\right].$

Note that t₀ = 0 in equations above. Equation 67 describes the time evolution of the fuel consumption during a long-time maneuver.

The solution of the two-point boundary value problem of going from an initial orbit O₀: (a₀, e₀, θ₀) to a final orbit O_f: (a_f, e_f, θ_f) defined by Eqs. 41-44 involves the determination of the initial value of the adjoint variables p'' _a0, p'' _φ0 and p'' _θ0, or equivalently, the determination of the auxiliary constants k₀, k₁ and k₂. Note that k₂ is written in terms of e₀, θ''₀ and k₁ through Eqs. 43 and 50. Accordingly, the solution of the two-point boundary value problem reduces to determine the other two constants.

Note that from the preceding section, Eq. 61, θ'' only depends on k₁. Accordingly, k₁ can be determined by solving Eq. 61, iteratively by means of Newton-Raphson method, for the given final conditions; that is, for θ''(t_f) = θ_f. Thus,

${\mathit{k}}_{1\mathit{n}+1}={\mathit{k}}_{1\mathit{n}}-\frac{{\mathit{\theta }}^{″}\left({\mathit{k}}_{1\mathit{n}}\right)-{\mathit{\theta }}_{\mathit{f}}}{{\left(\partial {\mathit{\theta }}^{″}/\partial {\mathit{k}}_1\right)}_{{\mathit{k}}_1={\mathit{k}}_{1\mathit{n}}}},$

with the partial derivative (∂θ''/∂k₁)_k1=k1n computed from Eq. 62. Note that the Newton-Raphson algorithm fails, if a conjugate point exists.

The constant k₀ is then determined as follows. Given k₁, the auxiliary variable ψ_f can be calculated from Eq. 55. Solving Eq. 41 for k₀, one finds:

(68) $\tan {\mathit{k}}_0=\frac{\sin \left(\sqrt{2}{\mathit{\psi }}_{\mathit{f}}\right)}{\sqrt{{\mathit{\alpha }}_{\mathit{f}}}-\cos \left(\sqrt{2}{\mathit{\psi }}_{\mathit{f}}\right)},$

Where ${\mathit{\alpha }}_{\mathit{f}}=\frac{{\mathit{a}}_{\mathit{f}}}{{\mathit{a}}_0}.$

Once the two-point boundary value problem has been solved, the fuel consumption variable J can be calculated from Eq. 67.

The initial value of the adjoint variables p''_a0, p''_φ0 and p_θ0 are then obtained as follows. From Eqs. 41 and 64, one finds:

(69) ${{\mathit{p}}^{″}}_{{\mathit{a}}_0}=\frac{{\mathit{v}}_0^2}{2{{\mathit{a}}^{″}}_0{\mathit{t}}_{\mathit{f}}}\left(\frac{{\mathit{u}}_{\mathit{f}}}{{\mathit{v}}_0}\right)\cos {\mathit{k}}_0,$

with u_f /ν₀ given by Eq. 66.

Now, solving Eqs. 48 and 49 for p''_φ0 and p''_θ0, it follows that

(70) ${{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2=\frac{8{\left({{\mathit{a}}^{″}}_0{{\mathit{p}}^{″}}_{{\mathit{a}}_0}\right)}^2}{{\cot }^2{\mathit{k}}_0\left(5{\sec }^2{\mathit{k}}_1-4\right)},$

(71) ${{\mathit{p}}^{″}}_{{\mathit{\phi }}_0}^2={{\mathit{p}}^{″}}_{{\mathit{\theta }}_0}^2\left({\sec }^2{\mathit{k}}_1-{\sec }^2{\mathit{\phi }}_0\right)$

The steps of an iterative algorithm for solving the two-point boundary value problem can be described as:

• For a given set of initial conditions (a₀, e₀) and final conditions (a_f, e_f), compute α_f, φ₀ = sin^-1e₀ and φ_f = sin^-1e_f.

• Guess a starting value for k₁.

• Compute τ₀ and τ_f from Eq. 43.

• Compute θ''(t_f) from Eqs. 44 and 50.

• If θ''(t_f) ≠ θ''_f, adjust the value of k₁ (through Newton-Raphson algorithm) and repeat steps 2 and 3 until |θ''(t_f) - θ''_f, | < δ, where δ is a prescribed small quantity.

• Compute ψ_f from Eq. 55.

• Compute (u_f/ν₀) from Eq. 66.

• Compute k₀ using Eq. 68.

• Compute successively p''_a0, p''_φ0 and p''_θ0 using the Eqs. 69-71.

The solution of the two-point boundary value problem for long duration transfers is used as starting guess for the solution of the complete problem as described later.

TRANSFERS BETWEEN CLOSE ORBITS

For transfers between close non-coplanar coaxial orbits, a simplified and complete first-order solution can be obtained through a linearization of the right-hand side of Eqs. 56-58, 59 and 60 about a reference orbit O - with semi-major axis a and eccentricity e. This simplified solution can be put in the form:

(72) $\Delta \mathit{x}={\mathit{Ap}}_0,$

where $\Delta \mathit{x}={\left[\Delta \mathit{\alpha }\Delta \mathit{e}\Delta \mathit{\theta }\right]}^T$ denotes the imposed changes on orbital elements (state variables), α = a/a, p_α = ap_a, p₀ denotes the 3 × 1 vector of initial values of the adjoint variables, and, A is a 3 × 3 symmetric matrix. In this simplified solution, the adjoint variables are constant and the matrix A is given by:

(73) $\mathit{A}=\left[\begin{array}{ccc}{\mathit{a}}_{\mathit{\alpha }\mathit{\alpha }}& {\mathit{a}}_{\mathit{\alpha }\mathit{e}}& {\mathit{a}}_{\mathit{\alpha }\mathit{\theta }}\\ {\mathit{a}}_{\mathit{e}\mathit{\alpha }}& {\mathit{a}}_{\mathit{ee}}& {\mathit{a}}_{\mathit{e}\mathit{\theta }}\\ {\mathit{a}}_{\mathit{\theta }\mathit{\alpha }}& {\mathit{a}}_{\mathit{\theta }\mathit{\alpha }}& {\mathit{a}}_{\mathit{\theta }\mathit{\theta }}\end{array}\right],$

where:

(74) ${\mathit{a}}_{\mathit{\alpha }\mathit{\alpha }}=4\sqrt{\frac{{\overline{\mathit{a}}}^5}{{\mathit{\mu }}^3}}{\left[\overline{\mathit{M}}+2\overline{\mathit{e}}\sin \overline{\mathit{E}}\right]}_{{\mathit{t}}_0}^{\mathit{t}}$

(75) ${\mathit{a}}_{\mathit{\alpha }\mathit{e}}={\mathit{a}}_{\mathit{e}\mathit{\alpha }}=4\sqrt{\frac{{\overline{\mathit{a}}}^5}{{\mathit{\mu }}_3}}{\left[\left(1-{\overline{\mathit{e}}}^2\right)\sin \overline{\mathit{E}}\right]}_{{\mathit{t}}_0}^{\mathit{t}}$

(76) ${\mathit{a}}_{\mathit{\alpha }\mathit{\theta }}={\mathit{a}}_{\mathit{\theta }\mathit{\alpha }}=0$

(77) ${\mathit{a}}_{\mathit{ee}}=\sqrt{\frac{{\overline{\mathit{a}}}^5}{{\mathit{\mu }}^3}}\left(1-{\overline{\mathit{e}}}^2\right){\left[\frac{5}{2}\overline{\mathit{M}}-\frac{5}{4}\overline{\mathit{e}}\sin \overline{\mathit{E}}+\frac{3}{4}\sin 2\overline{\mathit{E}}-\frac{1}{12}\overline{\mathit{e}}\sin 3\overline{\mathit{E}}\right]}_{{\mathit{t}}_0}^{\mathit{t}}$

(78) ${\mathit{a}}_{\mathit{e}\mathit{\theta }}={\mathit{a}}_{\mathit{\theta }\mathit{e}}=0$

(79) ${\mathit{a}}_{\mathit{\theta }\mathit{\theta }}=\sqrt{\frac{{\overline{\mathit{a}}}^5}{{\mathit{\mu }}^3}}\frac{1}{\left(1-{\mathit{e}}^2\right)}{\left[\frac{1+4{\overline{\mathit{e}}}^2}{2}\overline{\mathit{M}}+\left(-\frac{9}{4}\overline{\mathit{e}}+{\overline{\mathit{e}}}^3\right)\sin \overline{\mathit{E}}+\left(\frac{1}{4}+\frac{1}{2}{\overline{\mathit{e}}}^2\right)\sin 2\overline{\mathit{E}}-\frac{1}{12}\overline{\mathit{e}}\sin 3\overline{\mathit{E}}\right]}_{{\mathit{t}}_0}^{\mathit{t}}$

E - is the eccentric anomaly determined from Kepler's equation with the mean anomaly M = M₀ + n-(t - t₀) and t₀ is the initial time. Since the imposed variations on the orbital elements are linear in the adjoint variables, Eq. 72, the solution of the two-point boundary value problem is very simple and can be obtained by standard techniques.

For transfers between close orbits, the consumption variable can be written as:

(80) $\mathit{J}=\frac{1}{2}\left[{\mathit{a}}_{\mathit{\alpha }\mathit{\alpha }}{\mathit{p}}_{\mathit{\alpha }}^2+{\mathit{a}}_{\mathit{ee}}{\mathit{p}}_{\mathit{e}}^2+{\mathit{a}}_{\mathit{\theta }\mathit{\theta }}{\mathit{p}}_{\mathit{\theta }}^2\right]$

with a_αα, a_ee, and a_θθ given by Eqs. 74, 77 and 78, respectively.

SOLUTION OF THE TWO-POINT BOUNDARY VALUE PROBLEM

An iterative algorithm based on the complete first-order analytical solution is briefly described for solving the two-point boundary value problem of going from an initial orbit O₀:(a₀, e₀, θ₀) to a final orbit O_f: (a_f, e_f, θ_f) at the prescribed final time t_f.

For a given final time, Eqs. 56-58, 59 and 60 can be represented as follows:

${\mathit{y}}_{\mathit{i}}\left({\mathit{t}}_{\mathit{f}}\right)={\mathit{g}}_{\mathit{i}}\left({\mathit{t}}_{\mathit{f}},{\mathit{p}}_{{\mathit{a}}_0},{\mathit{p}}_{{\mathit{e}}_0},{\mathit{p}}_{{\mathit{\theta }}_0}\right),\mathit{i}=1,2,3,$

where y₁= a, y₂= e and y₃= θ. Note that p_a0, p_e0 and p_θ0 appear explicitly in the short periodic terms and also implicitly through a' (t), e' (t), θ'(t) and E'(t). Thus, functions g_i, i = 1, 2, 3, are nonlinear in these variables.

So, the two-point boundary value problem can be stated as: Find p_a0, p_e0 and p_θ0 such that y₁(t_f) = a_f, y₂(t_f) = e_f, y₃(t_f) = θ_f. This problem can be solved through a Newton-Raphson algorithm (Stoer and Bulirsch 2002), with the partial derivatives of functions g_i computed numerically by means of a procedure of centered differences. As previously mentioned, the starting guess for the iterative procedure is given by the solution of the two-point boundary value problem concerning transfers with long duration.

RESULTS

Two types of problems are considered: a direct problem, which corresponds to generate extremals trajectories for a given set of initial conditions for the state and adjoint variables, and an indirect problem concerning the solution of the two-point boundary value problem. In the first problem, a comparison between the complete, nonlinear and linear, first-order analytical solutions, derived in the preceding sections, and a numerical solution obtained by integrating the canonical system of differential equations governed by new Hamiltonian function H*, given by Eqs. 14 and 15, is discussed and the accuracy of the analytical approach is established. The second problem involves the comparison of the complete first-order analytical solution, including the short periodic terms, and the secular analytical solution (which describes the long duration transfers) in solving the two-point boundary value problem of going from an initial orbit O₀:(a₀, e₀, θ₀) to a final orbit O_f: (a_f, e_f, θ_f) at the prescribed final time t_f. A Runge-Kutta-Fehlberg method of orders 4 and 5, with step-size control, relative error tolerance of 10^-11, and absolute error tolerance of 10^-12, as described in Forsythe et al. (1977), has been used in the numerical integrations.

DIRECT PROBLEM

Figures 2-4 show the results for a comparison between three distinct solutions: the complete first-order analytical solution, the secular analytical solution, and the numerical solution. Two sets of initial conditions (state and adjoint variables) and transfer duration, defined in Table 1, are used in the comparison. Note that a transfer with moderate time of flight and a long duration transfer are computed. In Table 2, final values of state variables are shown. In these tables and in all figures, the state variables are expressed in canonical units, except the inclination of the orbital plane given in degrees. It should be noted that similar results can be obtained for maneuver with modification of the ascending node, since this maneuver is equivalent to the one with modification of the inclination of the orbital plane, taking into account the conditions described earlier (see paragraph after Eqs. 2-7). For simplicity, the results refer to maneuvers with modification of the orbital plane.

Figure 2
Comparison between secular, analytical and numerical time evolution of semi-major axis for maneuver 1 with t_f - t₀ = 25.0 and t_f - t₀ = 500.0.

Figure 3
Comparison between secular, analytical and numerical time evolution of eccentricity for maneuver 1 with t_f - t₀ = 25.0 and t_f - t₀ = 500.0.

Figure 4
Comparison between secular, analytical and numerical time evolution of inclination of orbital plane for maneuver 1 with t_f - t₀ = 25.0 and t_f - t₀ = 500.0.

Table 1
Set of initial conditions and transfer duration.

State and adjoint variables	Maneuver 1	Maneuver 2
a 0	1.00	1.00
e0	0.10	0.10
I0(degrees)	10.0	10.0
pa 0	4.90002 × 10-5	6.07832 × 10-4
pe 0	1.15518 × 10-5	1.92206 × 10-4
pI 0	1.28967 × 10-4	3.99208 × 10-4
	25	500

Table 2
Final state variables.

Maneuver	State variables	Analytical theory				Numerical solution
		Secular solution		Complete solution		Runge-Kutta-Fehlberg 4-5
	tf- t0	25	500	25	500	25	500
1	af	1.00491	1.10294	1.00489	1.10294	1.00489	1.10294
	ef	0.10071	0.11348	0.10066	0.11349	0.10066	0.11349
	If (degrees)	10.09730	12.05344	10.09684	12.05351	10.09684	12.05351
	Jf	2.3338 × 10-7	4.6678 × 10-6	2.3205 × 10-7	4.6678 × 10-6	2.3200 × 10-7	4.6680 × 10-6
2	af	1.06354	5.06122	1.06199	4.91725	1.06199	4.92809
	ef	0.11222	0.43929	0.10953	0.41415	0.10952	0.41691
	If (degrees)	10.31167	32.02676	10.31524	32.52541	10.31515	32.26476
	Jf	2.0662 × 10-5	4.1324 × 10-4	1.9995 × 10-5	4.1229 × 10-4	1.9986 × 10-5	4.1111 × 10-4

From the results presented in Figs. 2-4 and 5-7, and in Table 2, note that there exists an excellent agreement between the complete analytical solution and the numerical solution. The analytical solution has the same accuracy of the numerical solution for all state variables - semi-major axis, eccentricity, inclination of the orbital plane and consumption variable - considering moderate time of flight (t_f - t₀ = 25.0) and large time of flight (t_f - t₀ = 500.0). On the other hand, by analyzing the time evolution of the semi-major axis, eccentricity and inclination of the orbital plane represented in the plots of Figs. 2-4 and 5-7, as given by the complete analytical solution, the numerical solution and the secular solution, we can see that the amplitudes of the short periodic terms are small, but they can be significant for transfers with short or moderate duration. Note that the difference between the numerical solution and the secular solution decrease as the transfer duration increases. The good agreement between the numerical solution 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 previously described. This last remark is relevant, since the numerical algorithm to compute optimal low-thrust trajectories based on the analytical solution requires a simple Newton-Raphson method to solve the two-point boundary value problem, as described above.

Figure 5
Comparison between secular, analytical and numerical time evolution of semi-major axis for maneuver 2 with t_f - t₀ = 25.0 and t_f - t₀ = 500.0.

Figure 6
Comparison between secular, analytical and numerical time evolution of eccentricity for maneuver 2 with t_f - t₀ = 25.0 and t_f - t₀ = 500.0.

Figure 7
Comparison between secular, analytical and numerical time evolution of inclination of orbital plane for maneuver 2 with t_f - t₀ = 25.0 and t_f - t₀ = 500.0.

INDIRECT PROBLEM

In view of the previous results, one sees that the complete analytical first-order solution gives an accurate solution for the extremals trajectories concerning the transfers between coaxial non-coplanar orbits. So, in this section a study of some maneuvers is made using the analytical solution. The solution of the two-point boundary value problem of going from an initial orbit O₀: (a₀, e₀, θ₀) to a final orbit O_f: (a_f, e_f, θ_f) at the prescribed final time t_f is obtained considering the two versions of the analytical solution: the complete non-linear solution and the simplified linear solution for transfers between close orbits. This analysis allows us to discuss the applicability of the linear solution. Eight different maneuvers are considered, as described in Table 3. The initial orbit O₀ is the same for all maneuvers and it is defined by the following set of orbital elements: a₀ = 1.0 canonical unit, e₀ = 1.0 and I₀ = 1.0 degrees. The final value of the consumption variable J is taken as a comparison parameter for five transfer durations. Table 4 shows the results obtained through the two analytical solutions. J_Non-linear is the consumption variable computed by Eq. 60, and J_Linear is the consumption variable computed by Eq. 80. In both cases, variable J is expressed in canonical units. Figures 8-13 show the time evolution of the semi-major axis, the eccentricity, the inclination of the orbital plane and the consumption variable for maneuvers 1, 3 and 7. Two transfer durations are considered: t_f - t₀ = 25.0 and t_f - t₀ = 250.0 time units. All variables are expressed in canonical units, except the inclination of the orbital plane, given in degrees. It should be noted that there exists a good agreement between the two analytical solutions, linear and non-linear, for maneuver 1, which corresponds to a transfer between close orbits. On the other hand, for maneuver 3, which corresponds to a transfer with moderate amplitude, the results given by the linear solution are not good. Similar conclusion applies to maneuver 7, although the behavior for larger duration transfer is fair.

Table 3
Set of terminal orbits.

Maneuvers	Terminal orbits
	af (canonical units)	ef	If (degree)
1	1.05	0.105	11
2	1.50	0.45	12
3	2.00	0.25	15
4	1.10	0.15	11
5	1.60	0.25	11
6	1.50	0.30	12
7	1.50	0.20	15
8	1.50	0.25	15

Table 4
Consumption variable.

Maneuver	tf - t0	JLinear	JNon-linear	Maneuver	tf - t0	JLinear	JNon-linear
1	5	1.2257 × 10-04	1.2243 × 10-04	5	5	5.1744 × 10-03	6.0723 × 10-03
	10	5.7262 × 10-05	5.7226 × 10-05		10	2.2223 × 10-03	2.1756 × 10-03
	15	3.9004 × 10-05	3.9128 × 10-05		15	1.7213 × 10-03	2.0148 × 10-03
	20	2.8566 × 10-05	2.8568 × 10-05		20	1.1276 × 10-03	1.2010 × 10-03
	25	2.3348 × 10-05	2.3352 × 10-05		25	1.0282 × 10-03	1.1878 × 10-03
2	5	8.4019 × 10-03	8.7410 × 10-03	6	5	5.2947 × 10-03	5.8923 × 10-03
	10	3.2064 × 10-03	3.6288 × 10-03		10	2.0128 × 10-03	2.0770 × 10-03
	15	2.7985 × 10-03	2.9243 × 10-03		15	1.7472 × 10-03	1.9542 × 10-03
	20	1.7661 × 10-03	1.9309 × 10-03		20	1.0867 × 10-03	1.1831 × 10-03
	25	1.6413 × 10-03	1.7601 × 10-03		25	1.0112 × 10-03	9.9530 × 10-04
3	5	8.3527 × 10-03	1.0240 × 10-02	7	5	4.8712 × 10-03	5.2944 × 10-03
	10	5.3048 × 10-03	4.8075 × 10-03		10	2.1484 × 10-03	2.2786 × 10-03
	15	2.7761 × 10-03	3.4134 × 10-03		15	1.6252 × 10-03	1.7622 × 10-03
	20	2.4108 × 10-03	2.4108 × 10-03		20	1.1214 × 10-03	1.1969 × 10-03
	25	1.6888 × 10-03	2.0480 × 10-03		25	9.6155 × 10-04	9.6277 × 10-04
4	5	4.7730 × 10-04	4.9924 × 10-04	8	5	5.3652 × 10-03	5.8071 × 10-03
	10	1.8159 × 10-04	1.8234 × 10-04		10	2.2580 × 10-03	2.4075 × 10-03
	15	1.1375 × 10-04	1.1402 × 10-04		15	1.7871 × 10-03	1.9357 × 10-03
	20	9.3576 × 10-05	9.3779 × 10-05		20	1.1957 × 10-03	1.2827 × 10-03
	25	7.6080 × 10-05	7.6363 × 10-05		25	1.0520 × 10-03	1.0449 × 10-03

Figure 8
Time evolution of state variables for maneuver 1 with t_f - t₀ = 25.0.

Figure 9
Time evolution of state variables for maneuver 1 with t_f - t₀ = 250.0.

Figure 10
Time evolution of state variables for maneuver 3 with t_f - t₀ = 25.0.

Figure 11
Time evolution of state variables for maneuver 3 with t_f - t₀ = 250.0.

Figure 12
Time evolution of state variables for maneuver 7 with t_f - t₀ = 25.0.

Figure 13
Time evolution of state variables for maneuver 7 with t_f - t₀ = 250.0.

Now, consider the problem of solving the two-point boundary value problem for long-time transfers using the secular solution and the complete analytical solution, which includes the short periodic terms. In order to compare these solutions, maneuver 3 described in Table 3 is considered, with two transfer durations, t_f - t₀ = 125.0 and 500.0 time units. Table 5 shows the initial value of the adjoint variables obtained by using the secular solution and by using the complete solution. Note that there exist small differences between these initial values. As previously described, a small adjust of the initial value of the adjoint variables is made by a Newton-Raphson algorithm when the short periodic terms are included. Note that this adjust become smaller for a very large transfer duration. Moreover, Figs. 14-16 show that the secular solution of the TPBVP does not represent a mean solution of the complete solution of the same TPBVP.

Table 5
Initial value of the adjoint variables.

Maneuver		Secular solution			Complete solution
3	T	pa0	pe0	pI0	pa0	pe0	pI0
	125	0.001203	0.000370	0.000837	0.001215	0.000384	0.000801
	500	0.000300	0.000092	0.000209	0.000301	0.000093	0.000207

Figure 14
Secular solution and complete analytical solution of the TPBVP for maneuver 3 for both transfer durations - evolution of semi-major axis

Figure 15
Secular solution and complete analytical solution of the TPBVP for maneuver 3 for both transfer durations - evolution of eccentricity.

Figure 16
Secular solution and complete analytical solution of the TPBVP for maneuver 3 for both transfer durations - evolution of inclination of orbital plane.

Figure 17 represents the optimal trajectory for maneuver 3 with time of flight equal to 125.0 time units. Note that the motion of the vehicle resembles a spiral, departing from the initial orbit and arriving at the terminal orbit after twelve revolutions around the central body. If the maneuver represents a transfer between two orbits around the Earth with the semi-major axis of the initial orbit, a₀ = 1.0 distance unit, corresponding to 8000 km, then the maneuver lasts 39.35 hours and the magnitude of the average acceleration is equal to 1.59 cm/s². For such maneuver, the ratio between the average acceleration and the gravity acceleration on the ground is approximately 1.6 ' 10^-3, a typical value for low-thrust propulsion system. Similar conclusions are obtained for the other maneuvers discussed in this section.

Figure 17
Optimal trajectory for maneuver 3 with t_f - t₀ = 125.0.

CONCLUSION

An approximated first-order analytical solution for optimal time-fixed low-thrust limited power transfers (no rendezvous) between elliptic coaxial non-coplanar orbits in an inverse-square force field is determined through canonical transformation theory. The existence of conjugate point for long duration transfers is investigated through Jacobi condition. For transfers between close orbits a simplified solution is straightforwardly obtained by linearizing the new Hamiltonian and the generating function obtained through Hori method. For the direct problem of generating extremals trajectories, the analytical solution is compared to the numerical solution obtained by integrating the canonical system of differential equations describing the extremal trajectories for some sets of initial conditions, and the accuracy of the analytical approach is established. For the indirect problem, 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. A comparison between the linearized and nonlinear analytical solutions is made and it has been noticed that the linearized solution can provide good results for transfers between close orbits.

ACKNOWLEDGEMENTS

This research has been supported by CNPq under contract 304913/2013-8 and by FAPESP under contract 2012/21023-6.

REFERENCES

• Battin RH (1987) An introduction to the mathematics and methods of astrodynamics. New York: American Institute of Aeronautics and Astronautics, 796 p.
• Bliss GA (1946) Lectures on the calculus of variations. Chicago: University of Chicago Press.
• Bonnard B, Caillau JB, Dujol R (2006) Averaging and optimal control of elliptic Keplerian orbits with low propulsion. Systems & Control Letters 55(9):755-760. doi: 10.1016/j.sysconle.2006.03.004
  » https://doi.org/10.1016/j.sysconle.2006.03.004
• Camino O, Alonso M, Blake R, Milligan D, Bruin JD, Ricken S (2005) SMART-1: Europe's Lunar Mission Paving the Way for New Cost Effective Ground Operations (RCSGSO), Sixth International Symposium Reducing the Costs of Spacecraft Ground Systems and Operations (RCSGSO), European Space Agency ESA SP-601; Darmstadt, Germany.
• da Silva Fernandes S (2003) Notes on Hori method for canonical systems. Celestial Mechanics and Dynamical Astronomy 85(1):67-78. doi: 10.1023/a:1021711221662
  » https://doi.org/10.1023/a:1021711221662
• da Silva Fernandes S, das Chagas Carvalho F (2008) A first-order analytical theory for optimal low-thrust limited-power transfers between arbitrary elliptical coplanar orbits. Mathematical Problems in Engineering 2008. doi: 10.1155/2008/525930
  » https://doi.org/10.1155/2008/525930
• da Silva Fernandes S, das Chagas Carvalho F, de Moraes RV (2015) Optimal low-thrust transfers between coplanar orbits with small eccentricities. Computational and Applied Mathematics 1-14. doi: 10.1007/s40314-015-0249-9
  » https://doi.org/10.1007/s40314-015-0249-9
• Edelbaum TN (1964) Optimum low-thrust rendezvous and station keeping. AIAA Journal 2(7):1196-1201. doi: 10.2514/3.2521
  » https://doi.org/10.2514/3.2521
• Edelbaum TN (1965) Optimum power-limited orbit transfer in strong gravity fields. AIAA Journal 3(5):921-925. doi: 10.2514/3.3016
  » https://doi.org/10.2514/3.3016
• Edelbaum TN (1966) An asymptotic solution for optimum power limited orbit transfer. AIAA Journal 4(8):1491-1494. doi: 10.2514/3.3725
  » https://doi.org/10.2514/3.3725
• Edelbaum TN, Pines S (1970) Fifth and sixth integrals for optimum rocket trajectories in a central field. AIAA Journal 8(7):1201-1204. doi:10.2514/3.5872
  » https://doi.org/10.2514/3.5872
• Forsythe GE, Malcolm MA, Moler CB (1977) Computer methods for mathematical computations. Englewood Cliffs, New Jersey 07632. Prentice Hall, Inc., XI, 259 S. doi: 10.1002/zamm.19790590235
  » https://doi.org/10.1002/zamm.19790590235
• Geffroy S, Epenoy R (1997) Optimal low-thrust transfers with constraints-generalization of averaging techniques. Acta Astronautica 41(3):133-149. doi: 10.1016/s0094-5765(97)00208-7
  » https://doi.org/10.1016/s0094-5765(97)00208-7
• Gobetz FW (1964) Optimal variable-thrust transfer of a power-limited rocket between neighboring circular orbits. AIAA Journal 2(2):339-343. doi: 10.2514/3.2281
  » https://doi.org/10.2514/3.2281
• Gobetz FW (1965) A linear theory of optimum low-thrust rendezvous trajectories. Journal of the Astronautical Sciences 12:69.
• Haissig CM, Mease KD, Vinh NX (1993) Minimum-fuel, power-limited transfers between coplanar elliptical orbits. Acta Astronautica 29(1):1-15. doi: 10.1016/0094-5765(93)90064-4
  » https://doi.org/10.1016/0094-5765(93)90064-4
• Hori GI (1966) Theory of general perturbation with unspecified canonical variable. Publications of the Astronomical Society of Japan 18(4):287-296.
• Huang W (2012) Solving coplanar power-limited orbit transfer problem by primer vector approximation method. International Journal of Aerospace Engineering. doi: 10.1155/2012/480320
  » https://doi.org/10.1155/2012/480320
• Jamison BR, Coverstone V (2010) Analytical study of the primer vector and orbit transfer switching function. Journal of Guidance, Control, and Dynamics 33(1):235-245. doi: 10.2514/1.41126
  » https://doi.org/10.2514/1.41126
• Jiang F, Baoyin H, Li J (2012) Practical techniques for low-thrust trajectory optimization with homotopic approach. Journal of Guidance, Control, and Dynamics 35(1):245-258. doi: 10.2514/1.52476
  » https://doi.org/10.2514/1.52476
• Kiforenko BN (2005) Optimal low-thrust orbital transfers in a central gravity field. International Applied Mechanics 41(11):1211-1238. doi: 10.1007/s10778-006-0028-9
  » https://doi.org/10.1007/s10778-006-0028-9
• Lanczos C (1970) The variational principles of mechanics. 4th edition. Toronto: University of Toronto Press.
• Lie S (1888) Theorie der transformationsgruppen I. Leipzig: Teubner.
• Marec JP (1979) Optimal space trajectories. New York: Elsevier.
• Marec JP, Vinh NX (1980) Étude generale des transferts optimaux a poussee faible et puissance limitee entre orbites elliptiques quelconques. ONERA Publication 1980-1982.
• Marec JP, Vinh NX (1977) Optimal low-thrust, limited power transfers between arbitrary elliptical orbits. Acta Astronautica 4(5-6):511-540. doi: 10.1016/0094-5765(77)90106-0
  » https://doi.org/10.1016/0094-5765(77)90106-0
• Pines S (1964) Constants of the motion for optimum thrust trajectories in a central force field. AIAA Journal 2(11):2010-2014. doi: 10.2514/3.2717
  » https://doi.org/10.2514/3.2717
• Pontryagin LS, Boltyanskii VG, Gamkrelidze RV, Mishchenko EF (1962) The mathematical theory of optimal control processes. New York: John Wiley & Sons. 357 p.
• Quarta AA, Mengali G (2013) Trajectory approximation for low-performance electric sail with constant thrust angle. Journal of Guidance, Control, and Dynamics 36(3):884-887. doi: 10.2514/1.59076
  » https://doi.org/10.2514/1.59076
• Racca GD (2003) New challengers to trajectory design by the use of electric propulsion and other new means wandering in the solar system. Celestial Mechanics and Dynamical Astronomy 85(1):1-24. doi: 10.1023/a:1021787311087
  » https://doi.org/10.1023/a:1021787311087
• Racca GD, Marini A, Stagnaro L, van Dooren J, di Napoli L, Foing BH, Lumb R, Volp J, Brinkmann J, Grünagel R et al. (2002) SMART-1 mission description and developments status. Planetary and Space Science 50:1323-1336. doi:10.1016/S0032-0633(02)00123-X
  » https://doi.org/10.1016/S0032-0633(02)00123-X
• Rayman MD, Varghese P, Lehman DH, Livesay LL (2000) Results from the Deep Space 1 technology validation mission. Acta Astronautica 47(2-9):475-487. doi: 10.1016/s0094-5765(00)00087-4
  » https://doi.org/10.1016/s0094-5765(00)00087-4
• Stoer J, Bulirsch R (2002) Introduction to numerical analysis. 3rd edition. New York, USA: Springer.
• Sukhanov AA (2000) Optimization of low-thrust interplanetary transfers. Cosmic Research 38(6):584-587. doi: 10.1023/a:1026694802503.
  » https://doi.org/10.1023/a:1026694802503

Edited by

Section Editor: Antonio Bertachini

Publication Dates

• Publication in this collection
  2018

History

• Received
  31 Jan 2017
• Accepted
  01 June 2017