## Articulo

• Similares en SciELO

## versión impresa ISSN 0185-1101

### Rev. mex. astron. astrofis vol.54 no.1 Ciudad de México abr. 2018

Artículos

A numerical-analytical approach based on canonical transformations for computing optimal low-thrust transfers

1Departamento de Matemática, Instituto Tecnológico de Aeronáutica, Brazil.

2Instituto de Ciencia e Tecnologia, Universidade Federal de Sao Paulo, Brazil.

Abstract

A numerical-analytical procedure based on infinitesimal canonical transformations is developed for computing optimal time-fixed low-thrust limited power transfers (no rendezvous) between coplanar orbits with small eccentricities in an inverse-square force field. The optimization problem is formulated as a Mayer problem with a set of non-singular orbital elements as state variables. Second order terms in eccentricity are considered in the development of in the maximum Hamiltonian describing the optimal trajectories. The two-point boundary value problem of going from an initial orbit to a final orbit is solved by means of a two-stage Newton-Raphson algorithm which uses an infinitesimal canonical transformation. Numerical results are presented for some transfers between circular orbits with moderate radius ratio, including a preliminary analysis of Earth-Mars and Earth-Venus missions.

Key Words: celestial mechanics; methods: analytical; methods: numerical

Resumen

Se desarrolla un procedimiento numérico-analítico fundamentado en la transformación canónica infinitesimal para calcular transferencias de bajo empuje y potencia limitada con tiempo fijo entre órbitas coplanares de baja excentricidad en un campo de fuerza 1/r2. La optimización se formula como un problema del Mayer con un conjunto de elementos orbitales no singulares como variables de estado. Los términos de segundo orden en la excentricidad se consideran en el desarrollo del máximo Hamiltoniano que describe las trayectorias óptimas. El problema del valor l ́ımite de dos puntos al ir de una órbita inicial a una órbita final se resuelve mediante un algoritmo de Newton-Raphson de dos etapas que utiliza una transformación canónica infinitesimal. Se presentan resultados numéricos para transferencias entre órbitas circulares con una razón de radio moderada, y se incluye un análisis preliminar de las misiones Tierra-Marte y Tierra-Venus.

1. Introduction

The study presented in this paper was motivated by the use of low-thrust propulsion systems in new space missions. The two pioneer missions were NASA-JPL Deep Space 1 and ESA-SMART1. Deep Space 1 was the first interplanetary spacecraft to utilize Solar Electric Propulsion. It was launched on October 24, 1998, and its mission terminated on December 18, 2001, when its fuel supply was exhausted. SMART-1 was developed by ESA. 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 of the these pioneer missions can be found in Rayman et al. 2000; Racca et al. 2002; 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-4and 10-2) and show greatest benefits for high-energy planetary missions (Marec 1979; Racca 2003). Numerical and analytical solutions have been developed in the last fifty years 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 Carvalho, 2008; Jamison and Coverstone 2010; Jiang et al. 2012; Huang 2012; Quarta and Mengali 2013).

In previous works (da Silva Fernandes and das Chagas Carvalho 2008; da Silva Fernandes et al. 2016; das Chagas Carvalho et al. 2016) the authors presented complete first order analytical solutions based on canonical transformation theory for the problem of optimal time-fixed low-thrust limited-power transfers between arbitrary elliptic coplanar orbits, and for the problem of optimal time-fixed low-thrust limited-power transfers between coplanar orbits with small eccentricities. In the last study, non-singular orbital elements were used as state variables, and the maximum Hamiltonian and the optimal thrust acceleration were developed up to the first order in the eccentricity. It was observed that the analytical first order solution expressed in non-singular elements provides a good approximation when compared to a numerical solution obtained through a neighboring extremals algorithm considering large time of flight, moderate ratio of semi-major axis and arbitrary, very small, eccentricities. However, for a maneuver in which the radius of the final orbit is three times the radius of the initial orbit, the time behaviour of the eccentricity given by the analytical solution differs considerably from the behaviour given by the numerical solution. In that study, the numerical solution, used as exact solution, was obtained for the optimal transfer problem formulated as a Mayer problem with the radial distance and the components, radial and circumferential, of the velocity vector as state variables.

The main purpose of this paper is to develop a numerical-analytical procedure to solve the problem of optimal time fixed low-thrust limited power transfers (no rendezvous) between coplanar orbits with small eccentricities in an inverse-square force field, in which second order terms of the eccentricity are considered in the development of state equations such that solution can be applied to a larger class of maneuvers. The proposed procedure involves the development of a two-stage Newton-Raphson algorithm based on an infinitesimal canonical transformation. In the first stage of the algorithm, a two-point boundary value problem, defined by an average canonical system describing the secular behaviour of the optimal trajectories, is solved. The maximum Hamiltonian governing the average canonical system is derived by applying Hori (1966) method. So, an infinitesimal canonical transformation is built, and the short periodic terms are included in the solution by computing the Poisson brackets of the generating function with state variables. This technique provides an approximation to the numerical integration of the complete canonical system. However, when the short periodic terms are included in the solution of the two-point boundary value problem obtained in the first stage, small deviations from the prescribed final conditions arise, and the initial values of the adjoint variables computed in the first stage must be adjusted. The main advantage of this procedure is that the average canonical system is much simpler than the complete canonical system with periodic terms. Moreover, the average canonical system has a first integral, namely the new Hamiltonian, which is used to reduce the order of the system of differential equations and then to simplify the determination of the new values for the initial adjoint variables.

Numerical results are presented for the same transfers between circular coplanar orbits, including EarthVenus and Earth-Mars missions, discussed in the previous work (das Chagas Carvalho et al. 2016) and a comparison with the exact solution previously obtained is performed.

The paper is organized as follows. In § 2, the optimization problem for time-fixed low-thrust limited power transfer between coplanar orbits with small eccentricities is formulated as a Mayer problem of optimal control with a set of non-singular variables as state variables. In § 3, a numerical-analytical procedure to solve the two-point boundary value problem defined by the maximum Hamiltonian is described. In § 4, the Hori method is applied to eliminate the short periodic terms, and an infinitesimal canonical transformation is built. In § 5, the solution of the two-point boundary value problem of going from a initial orbit to a final orbit described by the average canonical system is presented. In § 6, numerical results are presented for some maneuvers, including Earth-Mars and Earth-Venus missions.

2. 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).

J=12t0tΓ2dt, (1)

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,

J=Pmax1m-1m0, (2)

where Pmax is the maximum power and m0 is the initial mass. The minimization of the final value of the fuel consumption, Jf , is equivalent to the maximization of mf .

Consider the motion of a space vehicle M, powered by a limited-power engine in an inverse-square force field. Since the terminal orbits are coplanar, the entire motion of the space vehicle is confined in a plane. So, at time t, the state of the space vehicle can be defined by the following set of non-singular orbital elements:

a,   h=ecosw,       k=esinw,     I=M+w, (3)

and the fuel consumption is defined by variable J. Here, a is the semi-major axis, e is the eccentricity, ω is the argument of pericenter and M is the mean anomaly.

The optimization problem can be formulated as a Mayer problem of optimal control as follows: It is proposed to transfer the space vehicle M from initial state (a0,h0,k0,l0,0) at time t0 = 0 to the final state (af , hf , kf , lf , Jf ) at time tf , such that the final consumption variable Jf is a minimum. The duration of the transfer tf − t0 is specified, and the final position of the vehicle is free, since only simple transfers (no rendezvous) are considered in this study.

The state equations, expressed up to the second order in eccentricity, are given by:

dhdt=1na2a-kcosl+14hsinl-118h2-138k2-1sinl-hcos2l+hsin2l-94hkcos3l+98h2-k2sin3lR+-32h-52h2-32k2-2cosl-hksinl+32hcosl+32ksin2l+32h2-k2cos3l+3hksin3lS (5)

dkdt=1nah138h2-118k2-1cosl+14hksinl-hcos2l-sin2l-98h2-k2cos3l-94hkcos3lR+-32k-hkcosl-32h2-52k2-2sinl-32kcos2l+32hsin2l-3hkcos3l+32h2-k2sin3lS (6)

dldt=n+1na-2-12h2+k2+32hcosl+32ksinl+12h2+k2cos2l+hksin2lR+-kcosl+hsinl+32hkcos2l+34h2+k2sin2lS (7)

dJdt=12(R2+S2) (8)

when n=μa3 is the mean motion, μ is the gravitational parameter, R and S are the radial and circumferential components of the thrust acceleration vector, respectively. Equations (4), (5), (6), (7) are obtained from the well-known Gauss equations (Battin 1987) using some expansions of the elliptic motion (Cayley 1861), up to second order in eccentricity, and the transformation of variables defined in Equation (3).

The performance index is defined by

IP=J(tf) (9)

For limited power system, it is assumed that there are no constraints on the thrust acceleration vector (Marec 1979).

According to the Pontryagin Maximum Principle (Pontryagin, et al, 1962), the adjoint variables pa, ph, pk, pl and pJ are introduced and the Hamiltonian function H(a, h, k, l, -, pa, ph, pk, pl, pJ , R, S) is formed:

H=paf1+phf2+pkf3+plf4+pJf5, (10)

where fi, i = 1, 2, 3, 4, 5 denotes the right-hand side of each differential equations in Equations (4), (5), (6), (7), (8), respectively.

The control variables, R and S, must be selected from the admissible controls such that the Hamiltonian function reaches its maximum along the trajectory. So, one finds

R*=-1napJ2a-kcosl+hsinl-2hkcos2l+h2+k2sin2lPa+-k-14hkcosl-118h2+138k2-1sinl-kcos2l-hsin2l-94hkcos3l+98h2+k2sin3lph+h+138h2+118k2-1cosl+14hksinl-hcos2l-ksin2l-98h2+k2cos3lpk+-2-12h2+k2+32hcosl+32ksinl+12h2+k2cos2l+hksin2lpi (11)

S*=-1napJ2a1-12h2+k2+hcosl+ksinl+h2+k2cos2l+2hksin2lPa+-32h-52h2+32k2-2cosl-hksinl+32hcos2l+32ksin2l+32h2+k2cos3l+3hksin3lph+-32-hkcosl-32h2+52k2-2sinl-32kcos2l+32hsin2l-3hkcos3l+32h2+k2sin3lpk+-kcosl+hsinl-32hkcos2l+34h2+k2sin2lpi (12)

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

From Equations (10), (11) and (12) one finds that the expression of the Hamiltonian, taking into account only terms up to the second order in eccentricty, is given by:

H*=H0+Hγ*, (13)

where

H0=npM

denotes the undisturbed Hamiltonian and Hγ * is the part related to the optimal thrust acceleration,

Hγ*=14n2a22+4ksinl+4hcosl+8hksin2l+4h2+k2cos2l4a2pa2+5-4k2-5h2+4ksinl-4hcosl-hksin2l+3-8h2-7k2cos2l+4hcos3l+4ksin3l+10hksin4l+5h2-5k2cos4lph2+5-5k2-4h2-4ksinl+4hcosl-hksin2l+7h2+8k2-3cos2l-4ksin3l-4hcos3l-10hksin4l+5h2-5k2cos4lpk2+-14h2+18k2-16sinl-4hkcos l+16sin2l-16kcos2l-18k2-h2sin3l-36hkcos 3lapapk+-4hksinl+16-18h2-14k2cosl+16ksin2l+16hcos2l+36hksin3l+18h2+k2cos3lapaph+-2hk-8hsinl-8kcosl+6-15h2+15k2sin2l+8hsin3l-8kcos 3l+10h2+10k2sin4l-20hkcos4lpkph+(32-48coslh-48sinlk-11cos2lh2+11k2cos2l+29h2+29k2-22hksin 2l)pi2+(8kcosl-8hsinl)apapi+-8sinl+kcos2l+7k-k2sin3l+3k2sinl-hsin2l-2hkcosl-2hkcos3l+h2sin3l+5h2sinlphpi+8cosl+ksin2l-5k2cos3l+k2cos3l-7h-hcos2l-2hksin3l+2hksin3l-h2cos3l-3h2coslpkpi (14)

It should be noted that in the development of Hγ *higher order terms in eccentricity are neglected, since such terms that arise from R*2 and S*2 are incomplete (they do not appear in R* and S*).

3. Numerical-analytical procedure to solve the two point-boundary value problem

In this section, a numerical-analytical procedure is described to solve the two-point boundary value problem of going from an initial orbit O0 at time t0 = 0 to a final orbit Of at a prescribed final time tf = T .

This two-point boundary value problem is defined by the canonical system of differential equations governed by the maximum Hamiltonian H* and the terminal constraints defining the initial orbit O0 and the final orbit Of ; that is,

dxdt=H*p,dpdt=H*x, (15)

with

at0=a0,                       ht0=h0,                     kt0=k0, (16)

atf=af,                       htf=hf,                     ktf=k0, (17)

In above equations, x denotes the state vector x = (a, h, k), p denotes de adjoint vector p = (pa, ph, pk), and are 3x1 matrices of partial derivatives of the maximum Hamiltonian with respect to the state vector and adjoint vector, respectively. For simple transfers, the final position of the space vehicle is free; so, variable l plays no special role in the two-point boundary value problem. Therefore, the differential equations for l and pl are not included in the set of differential equations above.

It should be noted that the explicit forms of the state and adjoint equations are relatively extensive due to the short periodic terms in the maximum Hamiltonian (see Eqn.(14)).

A two-stage Newton-Raphson algorithm is then proposed to solve the two-point boundary value problem described by Equations (15) (17). This algorithm has two distinct phases: in the first one, a two-point boundary value problem, defined by the average canonical system describing the secular behavior of the optimal trajectories, is solved. The maximum Hamiltonian governing the average canonical system is derived by applying the Hori (1966) method. So, an infinitesimal canonical transformation is built, and the short periodic terms are included in the solution by computing the Poisson brackets of the generating function with the state variables. This technique provides an approximation, to first order in a small parameter closely related to the magnitude of the optimal thrust acceleration, to the time behavior of the state variables considering the complete maximum Hamiltonian. In other words, it provides an approximation to the solution obtained by numerical integration of the complete canonical system. However, when the short periodic terms are included in the solution of the two-point boundary value problem obtained in the first stage, small deviations from the prescribed final conditions arise. Accordingly, the initial values of the adjoint variables computed in the first stage must be adjusted by a second Newton Raphson algorithm. The advantage of this procedure is that the average canonical system is much simpler than the complete system and it has a first integral used to reduce the order of the system of differential equations, making easier the numerical integration and the determination of the initial value for the adjoint variables.

In the next section, the Hori method is applied to obtain the infinitesimal canonical transformation and the average canonical system. In § 5, the solution of the two-point boundary value problem described by the average canonical system is briefly defined.

4. Infinitesimal canonical transformation

In order to eliminate the short periodic terms from the maximum Hamiltonian H*, an infinitesimal canonical transformation is built by means of the Hori (1966) method. It is assumed that the undisturbed Hamiltonian H0 is of zero order and the remaining part H* is of first order in a small parameter associated to the magnitude of the optimal thrust acceleration (da Silva Fernandes and das Chagas Carvalho, 2008; da Silva Fernandes et al, 2016; das Chagas Carvalho et al, 2016).

Consider an infinitesimal canonical transformation,

(a,h,k,l,pa,ph,pk,pl)(a´,h´,k´,l´,p´a,p´h,p´k,p´l)

The new variables are denoted by the symbol “´”.

Following the algorithm of the Hori method, one finds from the zero-th order equation of the algorithm that the new undisturbed Hamiltonian is given by

F0=H0, (18)

with F0 expressed in terms of the new canonical variables.

Now, consider the canonical system described by F0:

The general solution of this system is very simple and it is given by

a´=a0´,P´a=P0´+3n´(t-t0)2a´P´I,h´=h0´,P´h=Ph0´,k´=k0´,P´k=Pk0´,l´=l0´+n´(t-t0),P´l=Pl0´, (19)

The subscript ’0’ denotes the constants of integration.

Introducing the general solution defined above into the equation of order 1 of the algorithm of the Hori method results in

S1t=Hγ-F1. (20)

Following the integration theory for the Hori method proposed by da Silva Fernandes (2003), one finds that the new (average) Hamiltonian function F1 and the generating function S1 are given, respectively, by:

F1=12n2a24a2pa2+52Ph2+Pk2-hph+kpk2-2hph-kpk2 (21)

and

S1=12n2a242h2-2k2sin2l-4kcosl-4hkcos2l+4hsinla2pa2+-54h2+54k2sin4l+-4h2-72k2+32sin2l-4kcosl-4hsinl+12hkcos2l-52hkcos4l+43hsin3l-43kcos3lph2+72h2+4k2-32sin2l+-54h2+54k2sin4l+52hkcos4l+12hkcos2l+4hsinl+4kcosl-43kcos3lpk2+2-9h2-7k2+8sinl+3h2+3k2sin3l+4hsin2l 2hkcosl-4kcos2l-6hkcos3lapaph+2-3h2+3k2cos3l+7h2+9k2-8cosl-4hcos2l-2hksinl-6hksin3l-4ksin2lapapk+152h2-3+152k2cos2l+52k2-52h2cos4l-83ksinl+8hcosl-5hksin4l-8ksinl-83hcos3lphpk (22)

Terms factored by pl i have been omitted in equations above, since only transfers (no rendez-vous) are considered. Note that pl i is a first integral of the average canonical system; that is, pl i is constant and its value is computed from the transversality condition which gives pl(tf) = 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.

By computing the Poisson bracket of the generating function with the state variables a, h, k -, one finds that the infinitesimal canonical transformation built by means of the Hori method is then defined by the following equations:

a=a´+14(n´3a´2)8(2h´2-2k´2)sin 2l´-4k´cos l´-4h´k´cos 2l´+4h´sin l´a´2p´a+2a´(-9h´2-7k´2+8)sin l´+(3h´2-3k´2)sin 3l´+4h´sin 2l´+2h´k´cos l´-6h´k´cos 3l´-4k´cos 2l´p´h+2a´(-3h´2-3k´2)cos 3l´+(7h´2+9k´2-8)cos l´-4h´ cos 2l´-2h´k´sin l´-6h´k´sin 3l´-4k´cos 2l´P´k (23)

h=h´+14(n´3a´2)2(-54k´2+54h´2)sin 4l´+-4h´2-72k´2+32sin2l´-4k´cos l´-4h´sin l´+12h´k´ cos 2l´-52h´k´ cos 4l´+43h´sin3l´P´h+2a´-9h´2-7k´2+8sinl´+(3h´2-3k´2)sin 3l´+4h´ sin 2l´+2h´k´ cos l´-6h´k´ cos 3l´-4k´cos 2l´P´a+153h´2-3+152k´2cos2l´+(-52k´2-52h´2)cos4l´-83k´sin3l´+8h´cosl´-5h´k´sin4l´-8k´sinl´-83h´cos3l´P´k (24)

k=k´+14(n´3a´2)2(72h´2+4k´2-32)sin 2l´+-54h´2-54k´2sin4l´+52h´k´ cos 4l´+12h´k´ cos 2l´+4h´ sin l´+4k´cos l´-43h´ sin 3l´P´k+2a´-3h´2+3k´2cos 3l´+(7h´2+9k´2-8)cos l´-4h´ cos 2l´+2h´k´ sin l´-6h´k´ sin 3l´-4k´sin 2l´P´a+152h´2-3+152k´2cos2l´+(52k´2-52h´2)cos4l´-83k´sin3l´+8h´cosl´-5h´k´sin4l´-8k´sinl´-83h´cos3l´P´h (25)

Equations (23), (24) and (25) can be applied to get an approximation of the time evolution of the state variables if the solution of the average canonical systems governed by the new Hamiltonian F1 is known. In other words, they provide an approximation to the numerical solution obtained by integrating the complete canonical system governed by the maximum Hamiltonian H*. Note that l´ is given by numerical integration of the differential equation

dl´dt=n´+74a´μk´p´h-h´p´k. (26)

The second term of right-hand side of Equation (26) is obtained from the terms factored by pl in the new Hamiltonian, that dissapear when one takes pl = 0 for the case of simple transfers. Recall that pl is a constant, whose value obtained from the transversality condition for simple transfers is zero.

Equations (23) (25) can be put in the following compact form

a=a´+δa,h=h´+δh,k=k´+δk, (27)

where δa, δh and δk denotes the periodic terms in the right-hand side of Equations (23), (24) and (25). Applying the initial conditions one finds that the time evolution of the state variables is written in a compact form as

at=a´t+δat-δat0,ht=h´t+δht-δht0,kt=k´t+δkt-δkt0, (28)

where δa(t), δh(t) and δk(t) are calculated at time t, with a´,h´, k´, pa´, ph´, pk´ obtained by numerical integration of the average canonical system from the initial conditions a0, h0, k0, pa0 , ph0 and pk0 . δa(t0), δh(t0), δk(t0) are calculated in terms of the initial conditions at t = t0. Accordingly, the solution with periodic terms and the average solution satisfy the same initial conditions.

Note that if the two-point boundary value problem governed by the average canonical system is solved by means of a Newton-Raphson algorithm, that is, pa0 , ph0 and pk0 are determined such that the final conditions are satisfied, then the periodic terms expressed by δa(t) − δa(t0), δh(t) − δh(t0) and δk(t) − δk(t0) provide small deviations from the final conditions, and pa0 , ph0 , pk0 must be adjusted in order to satisfy these conditions. This adjustment is made by another Newton-Raphson algorithm in the second stage of the proposed algorithm.

5. The average two-point boundary value problem

In what follows is described the first stage of the proposed algorithm, in which the two-point boundary value problem defined by the average canonical system is solved. The Newton-Raphson algorithm developed in this first stage uses a first integral generated by the Hori method, namely the new Hamiltonian, to reduce the order of the canonical system and to simplify the determination of new values for the initial adjoint variables.

The two-point boundary value problem defined by the average canonical system governed by the new Hamiltonian F1 is then given by the following set of differential equations

with the boundary conditions given by

a0= a0,     atf=af,h0=h0,       htf=hf,k0=k0,       ktf=kf, (30)

As before, the symbols denoting the new canonical variables are omitted.

According to the Hori method, the new Hamiltonian F1 is a first integral of the canonical system defined above and can be used to reduce the order of the system. From the differential equations for a and pa one finds = −E.

d(apa)dt=-E. (31)

Thus,

atpat=a0pa0-ET. (32)

The constant E denotes the value of new Hamiltonian, that is, F1 = E.

Finally, using the above result, the differential equation for a can be integrated. After simplifications, one finds

at=a01+2a0μEt2-2Bt, (33)

where B = a0pa0 .

Equations (31) and (33) can be used to simplify the algorithm of the Newton-Raphson method, since only four differential equations have to be integrated. Note that three unknowns - pa0 , ph0 , pk0 - are iterated at each step of the algorithm.

We note that the average canonical system has two other first integrals:

kph-hpk=C1, (34)

Ph2+Pk2-kph+hpk2=C22. (35)

But, they do not allow further simplifications to the algorithm.

6. Results

In this section, the numerical-analytical procedure based on canonical transformation proposed to compute optimal time-fixed low-thrust limited power transfers (no rendezvous) between coplanar orbits with small eccentricities in an inverse-square force field is applied in a preliminary analysis of some transfers discussed in the previous work (das Chagas Carvalho et al. 2016).

As noted in the previous work, the complete first order analytical solution, with terms of first order in eccentricity, does not provide accurate results for the time behaviour of the eccentricity along the optimal trajectory, considering transfers between circular coplanar orbits with moderate amplitude (radius of the final orbit three times the radius of the initial orbit). Accordingly, the present analysis considers some values of the radius ratio ρ=rfr0(where r0 is the radius of the initial circular orbit and rf is the radius of the final circular orbit), and two values of the time of flight T = tf − t0. Table 1 shows values of the consumption variable J for five values of ρ and two values of T . In this table, J1 denotes the consumption variable computed by solving the transfer problem formulated as a Mayer problem of optimal, with the radial distance and, the radial and circumferential components of the velocity vector as state variables, and with the two-point boundary value problem solved through a neighbouring extremals algorithm based on the state transformation matrix (see das Chagas Carvalho et al. 2016 for more details). This approach of the transfer problem will be referred to as Method 1. J2 denotes the consumption variable computed by the proposed algorithm, which will be referred to as Method 2. Note that ρ = 0.7270 DU, ρ = 1.5236 DU and ρ = 2.5000 DU correspond to interplanetary missions f/or transfers from Earth to Venus, from Earth to Mars and from Earth to Asteroid Belt, respectively. In the preliminary interplanetary missions analysis, the following assumptions are considered:

Table 1 Consumption variable J

ρ=rfr0 Tf-t0 J1 J2
0.7270 25.0 5.9852 × 10−4 6.0010 × 10−4
125.0 1.1949 × 10−4 1.1950 × 10−4
1.5236 25.0 7.2468 × 10−4 7.3362 × 10−4
125.0 1.4421 × 10−4 1.4428 × 10−4
2.0000 100.0 4.2976 × 10−4 4.3134 × 10−4
200.0 2.1462 × 10−4 2.1485 × 10−4
2.5000 100.0 6.7826 × 10−4 6.8525 × 10−4
200.0 3.3811 × 10−4 3.3894 × 10−4
3.0000 100.0 9.0260 × 10−4 9.2706 × 10−4
200.0 4.4776 × 10−4 4.4505 × 10−4

1. The orbits of the planets are circular;

2. The orbits of the planets lie in the plane of the ecliptic;

3. The flight of the space vehicle takes place in the plane of the ecliptic;

4. Only the heliocentric phase is considered; that is, the attraction of planets on the spacecraft is neglected.

All results in the Table 1 are expressed in canonical units. The starting guess for the adjoint variables is given by the first order solution in eccentricity developed in previous works (da Silva Fernandes et al. 2016; das Chagas Carvalho et al. 2016). Note that the consumption variable J2 is computed as described in previous works; that is,

J2=Etf-t0+S1, (36)

where ∆S1 = S1(tf ) − S1(t0), with, S1 given by Equation (22).

It should be noted that the numerical results computed through the proposed algorithm (Method 2) for the consumption variable have almost the same accuracy as those of the numerical approach based on the neighbouring extremals algorithm (Method 1). This good agreement between the results is closely related to the excellent accuracy obtained by the proposed algorithm (Method 2) for the time behavior of the semi-major axis, which represents the main parameter in transfers between coplanar circular orbits (the fuel is mainly spent to satisfy the imposed change in the semi-major axis).

The time evolutions of semi-major axis and eccentricity are presented in Figures 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10; as well as the difference between the solutions obtained by the methods. From Figures 1 to 10 and Table 1, note that the excellent agreement between the results for all variables semi-major axis, eccentricity and consumption variable for all values of ρ and T presented in Table 1. The solution considering second order terms in eccentricity in the development of the state equations obtained by means of the proposed algorithm (Method 2) has almost the same accuracy as Method 1: in the worst case (ρ = 3.0000 DU and tf − t0 = 100.0 TU) the difference for the semimajor axis does not exceed 10-2DU, and the difference for the excentricity does not exceed 8 × 10-3.

In Figures 1-10 the contribuition of the periodic terms in the time behavior of the orbital elements, mainly, of the eccentricity should be noted. The analytical-numerical solution given by Equations (23), (24), (25) provides a good description for the solution of the two-point boundary value problem governed by the maximum Hamiltonian defined by Equations (13) and (14).

The optimal trajectories with amplitudes ρ = 0.7270 DU and ρ = 1.5236 DU and time of flight tf − t0 = 25.0 TU are depicted in Figure 11, and the optimal trajectories with amplitudes ρ = 2.0000 DU and ρ = 2.5000 DU and time of flight tf − t0 = 100.0 TU are depicted in Figure 12. Finally, the optimal trajectory with largest amplitude ρ = 3.0000 DU and time of flight tf − t0 = 100.0 TU is shown in Figure 13. As mentioned before, the optimal trajectories computed through the proposed algorithm (Method 2) give a good description for the solution of the two-point boundary value problem of going from an initial orbit to a final orbit when compared to the trajectories computed by means of the numerical approach based on the neighboring extremals algorithm (Method 1).

7. Conclusion

In this paper, a numerical-analytical procedure for computing optimal time-fixed low-thrust limited power transfers between coplanar orbits with small eccentricities in an inverse-square force field is presented. Second order terms in eccentricity are considered in the development of the state equations. The proposed procedure is based on an infinitesimal canonical transformation built by applying Hori method. A two-stage NewtonRaphson algorithm is described to solve the two-point boundary value problem of going from an initial orbit to a final orbit. The proposed algorithm is applied in a preliminary analysis of some interplanetary missions discussed in previous work. Numerical results show the very good agreement between the results provided by the proposed algorithm and the numerical results previously obtained through a different approach of the transfer problem. The time behavior of the eccentricity along the optimal trajectory shows excellent agreement with previous results, including maneuvers with larger radius ratio and time of flight.

Finally, we note that the procedure described in this paper can be extended to more general low-thrust limited power transfers.

References

Battin, R. H. 1987, An introduction to the mathematics and methods of astrodynamics, (2nd ed., New York, NY: AIAA) [ Links ]

Bonnard, B., Caillau, J. B., & Dujol, R. 2006, Averaging and optimal control of elliptic Keplerian orbits with low propulsion, Systems and control letters, 55, 755 [ Links ]

Camino, O., Alonso, M., Blake, R. et al. 2005, ESA SP601, Reducing the Costs of Spacecraft Ground Systems and Operations, (Darmstadt, Germany, ESA) 43 [ Links ]

Cayley, A. 1861, MmRAS, 29, 191 [ Links ]

da Silva Fernandes, S. 2003, CeMDA, 85, 67 [ Links ]

da Silva Fernandes, S. & das Chagas Carvalho, F. 2008, Mathematical Problems in Engineering, ID 525930, dx.doi.org/10.1155/2008/525930 [ Links ]

da Silva Fernandes, S., das Chagas Carvalho, F., & de Moraes, R. V. 2016, Optimal low-thrust transfers between coplanar orbits with small eccentricities, Computational and Applied Mathematics, 35, 803 [ Links ]

das Chagas Carvalho, F., da Silva Fernandes, S., & de Moraes, R. V. 2016, A numerical study for optimal lowthrust limited power transfers between coplanar orbits with small eccentricities, Computational and Applied Mathematics , 35, 907 [ Links ]

Edelbaum, T. N. 1964, AIAAJ, 2, 1196 . [ Links ]

______. 1965, AIAAJ, 3, 921 [ Links ]

__________. 1966, AIAAJ, 4, 1491 [ Links ]

Geffroy, S. & Epenoy, R. 1997, AcAau, 41, 133 [ Links ]

Gobetz, F. W. 1964, AIAAJ , 2, 339 . [ Links ]

______1965, JAnSC, 12, 69 [ Links ]

Haissig, C. M., Mease, K. D., & Vinh, N. X. 1993, AcAau , 29, 1 [ Links ]

Hori, G. I. 1966, PASJ, 18, 287 [ Links ]

Huang, W. 2012, International Journal of Aerospace Engineering, http://dx.doi.org/10.1155/2012/480320 [ Links ]

Jamison, B. R. & Coverstone, V. 2010, JGCD, 33, 235 [ Links ]

Jiang, F., Baoyin, H., & Li, J. 2012, JGCD , 35, 245 [ Links ]

Kiforenko, B. N. 2005, IAM, 41, 1211 [ Links ]

Marec, J. P. 1979, Optimal space trajectories, (Amsterdam, Netherland: Elsevier) [ Links ]

Marec, J. P. & Vinh, N. X. 1977, AcAau , 4, 511 [ Links ]

__________. 1980, E´tude generale des transferts optimaux a poussee faible et puissance limitee entre orbites elliptiques quelconques, ONERA Publication 1980 [ Links ]

Pontryagin, L. S., Boltyanskii, V. G., Gamkrelidze, R. V., & Mishchenko, E. F. 1962, The mathematical theory of optimal control processes, (New York, NY: John Wiley and Sons) [ Links ]

Quarta, A. A. & Mengali, G. 2013, JGCD , 36, 884 [ Links ]

Racca, G. D. 2003, CeMDA , 85, 1 [ Links ]

Racca, G. D, Marini, A., Stagnaro, L., et al. 2002, P&SS, 50, 1323 [ Links ]

Rayman, M. D, Varghese, P., Lehman, D. H., et al. 2000, AcAau , 47, 475 [ Links ]

Sukhanov, A. A. 2000, CosRe, 38, 584 [ Links ]

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

Received: September 06, 2017; Accepted: November 03, 2017

Sandro da Silva Fernandes and Francisco das Chagas Carvalho:

Departamento de Matemática, Instituto Tecnológico de Aeronáutica, Brazil.

João Victor Bateli Romão:

Instituto de Ciencia e Tecnologia, Universidade Federal de São Paulo, Brazil.