SciELO - Scientific Electronic Library Online

 
vol.55 número2Enhanced mass loss rates in red supergiants and their impact on the circumstellar mediumThe environmental dependence of galaxy age and stellar mass in the redshift region 0.6 ≤ z ≤ 0.75 índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

Links relacionados

  • No hay artículos similaresSimilares en SciELO

Compartir


Revista mexicana de astronomía y astrofísica

versión impresa ISSN 0185-1101

Rev. mex. astron. astrofis vol.55 no.2 Ciudad de México oct. 2019  Epub 25-Nov-2020

https://doi.org/10.22201/ia.01851101p.2019.55.02.05 

Articles

Minimizing fuel consumption in orbit transfers using universal variable and particle swarm optimization

L. M. Echeverry1  2 

Y. Villanueva1  3 

1School of Exact Sciences and Engineering, Sergio Arboleda University, Bogotá, Colombia.

2Department of Mathematics, University of the Andes, Bogotá, Colombia.

3Institute of Mathematics and Statistics, Federal University of Goiás, Goiania, Brazil.


ABSTRACT

Minimizing fuel consumption in space travels is becoming increasingly important for spatial development. In the present paper, the fuel consumption in orbit transfers (without gravitational assistance) is minimized, where a spacecraft performs a change from an orbit around the Earth to another one around a different celestial body. Two methods are presented: one of immediate transfer and another with wait time. Minimizing is done by solving a nonlinear system, obtained by applying Lagrange multipliers to the equation modelling the keplerian system, and using the seeds coming from the particle swarm algorithm to execute the Newton’s method. Numerical simulations with real values were made to compare these methods with the Hohmann transfer and data from the specialized literature.

Key Words: celestial mechanics; methods; numerical; planets and satellites; fundamental parameters; space vehicles

RESUMEN

La minimización del gasto de combustible en los viajes espaciales es cada día más importante para el desarrollo espacial. En el presente trabajo se minimiza el gasto de combustible en transferencias de orbita (sin asistencia gravitacional), donde se ejecuta un cambio de órbita de una nave alrededor de la Tierra a una órbita alrededor de otro cuerpo celeste. Se presentan dos métodos, uno de transferencia inmediata y otro con tiempo de espera. La minimización se hace resolviendo un sistema no lineal que aparece después de aplicar multiplicadores de Lagrange a las ecuaciones que modelan el sistema kepleriano, usando las semillas que vienen del algoritmo de enjambres de partículas para ejecutar el método de Newton. Se hicieron simulaciones numéricas con valores reales para comparar estos métodos con la transferencia de Hohmann y los datos que aparecen en la literatura especializada.

1. INTRODUCTION

Recently, authors such as (Gang et al. 2014; Shan et al. 2014; Zotes et al. 2012) have given important contributions to the study of spatial trajectory optimization. The latter have considered a geometric method with results far from those presented in this article, by having a large flight time with a v similar to ours.

Thus, to minimize fuel consumption in the acceleration v [as in (Sharaf & Saad 2016)], the present article considers the methods proposed in Leeghim (2013), which are variations of the Lambert problem, with neither direction of motion in a plane nor the time t of the transfer at the beginning given; both are calculated. Another difference of our approach is the obtainment of the wait time t1 (Vallado 1997), where the interceptor does not leave its initial orbit until the relative positions of the bodies are convenient. This is presented in § 2.

The problem becomes a constrained optimization system. Its direct solution leads to several systems of nonlinear equations (9 equations as in § 3, 11 equations as in § 4). The solution method, presented in § 5 (Newton’s method) needs an appropriate starting point (seed). This seed is usually found from an exhaustive model that takes a long time to locate a suitable value (grids). Therefore, an heuristic technique, the particle swarm optimization presented in § 5.1, is used to approximate the seed in a reasonable time. Such methods were applied to trips with real data in the solar system and are presented in § 6. The results of the first case (immediate transfer) were compared with information released from missions of several space agencies in Villanueva (2013) and here the two problems are compared with the data shown in Kemble (2006), who deals with the Lambert problem directly from a geometric point of view.

2. PRELIMINARIES

The two-body problem (Keplerian) described in Bate (1971) was considered. The equation

r¨t=-µrt3rt (1)

models the dynamic of the system, with µ = GM . The position of the particle on an ellipse is r(t) = [x1(t) x 2(t)], such that

x1t=acosEt-ϵ,             x22t=a1 - ϵ2 sinEt,

where a is the semi-major axis, ǫ the eccentricity, E the eccentric anomaly. Using the fact that the norm of the angular momentum of the particle is constant we get

E =µa3/21-ϵcosE,

and therefore Kepler’s equation is obtained:

Et-ϵsinEt=µa3/2t-t0=M, (2)

where M is the mean anomaly (Bate 1971). To find a solution E for Kepler’s equation (2) at any time normally Newton’s method is used, but it does not converge or converges too slowly when ǫ ≈ 1 (Elipe et al. 2017). From Danby (1962, p. 168) we use Taylor’s expansion in (2) to derive the equations of the dynamics related to the problem for any type of trajectory (classical curves). Then

µt-t0&=a3/2E- ϵ  E- E33!+E55!-E77!+· · ·=a3/2&1-ϵE+ϵ+aE33!-13!aE55! +1a2  aE77!-  ···.

With x = aE, we obtain

µt-t0&=a1ϵx+ϵx33!- 1-x55!+1-x77!- · · ·.

This equation accepts any value of 𝜖 and ȧ ≠ 0. Therefore, using the fact r = a(1 − ȧ 𝜖 cos E), (Bate 1971, p. 187), and with the expression of Ė, an universal variable, x ∈ ℝ is defined as := µr .

So, the following expressions are obtained for t and r (Bate 1971) using equation (1), integrating the universal variable as µdt = rdx:

µt = a x-asinxa+r0 ,v0µa  1-cosxa +r0 asinxa, (3)

r = a+a r0 ,v0µasinxa+r0a-1cosxa.

Since the initial position and velocity vectors r0 and v0 are linearly independent, and having r, and r = r in the same plane, the vectors of position and velocity are expressed in terms of the initial vectors and the universal variable (Bate 1971):

rt=fr0+gv0v,vt=f˙r0+g˙v0. (4)

3. PROBLEM STATEMENT

The target travels in its own orbit, and the interceptor orbits around the Earth. To calculate the orbit to be taken by the interceptor to fly by the target the initial position and initial velocity of the interceptor r0, v0 and the target R0, V0 are required. We look for the minimum change of velocity Δv0, vtransfer = v0 + Δv0, allowing overflight to be achieved. See Figure 1.

Fig. 1 Description of the first problem. 

Using equation (3), let A and X be the semi-major axis and the universal variable, respectively, associated with R, and a and x the semi-major axis and the universal variable, respectively, associated with r (in the orbit transfer). The following functions modelling the movement of the two bodies in the solar system are defined:

η1=AX-A sin XA+ R0 V0µA1 - cos XA+R0A sinXA-µt,

η2=ax-a sin xa+ r0 v0+Δv0µa1 - cos xa+r0a sinxa-µt.

We call ηs: = (η1 η 2). Let R and rbe the positions of the target and the interceptor, respectively, shown in equation (4) and let the functional J = 1/2 Δv0 T Δv0. Then we have the following optimization problem:

  • Minimize J(v 0 ).

  • Restricted to ηs(X, x, t, Δv0) = 0 and (Rr)(X, x, t, Δv0) = 0.

To solve the problem, Lagrange multipliers are used in the following functional:

Hs=JΔv0+λTηs+ϕTR-r,

where λ ∈ ℝ2 and ϕ ∈ ℝ 3.

After setting ∇Hs = 0 as Leeghim (2013), we obtain the following system of nonlinear equations:

fs =ηs=η1,η2=0,R -r= 0,ϕTRX-rx+rµRt-rt= 0,J  Δv0+ 1r ϕT r x η2 Δv0 - ϕTr Δv0 = 0. (5)

4. WAIT TIME

From Vallado (1997, p. 318) the wait time for planar and circular orbits is given by

t1 = θ - θi + 2kπW-w

where θi and θ are the initial and final angles (after t1) between R and r , and w and W are the angular speeds of the interceptor and the target, respectively. In this paper we calculate t1 (the wait time), which is one of the results of the problem to be minimized for any kind of transfer orbit (in space) and maintains the same functional to be minimized. See Figure 2.

Let ηc = (η1 η2 η3) such that:

η1=A X-A sin XA+R0V0µ A 1 - cos XA+R0 A sin- xA- µt,

η2 =a0  x1-a0sinx1a0+r0 v0µ a01-cosx1a0+r0a0sinx1a0-µt1,

η3 =a  x-asinxa+r1 v1+Δ v1µ a01-cosxa+r1a0sinxa0-µt-t1,

Fig. 2 Description of the second problem. 

where η1 is the equation of motion of the target, η2 the equation of motion for the interceptor in its initial orbit during the wait time, and η3 is the equation of motion of the interceptor after the wait time. Then we have a problem similar to the previous one:

  • Minimize J(Δv1).

  • Restricted to ηc(X, x1, x, t, t1, Δv1) = 0 and (Rr )(X, x1, x, t, t1, Δv1) = 0.

Where J =1/2 Δv1 T Δv1. Therefore, the new functional is:

Hc=JΔv1+λTηc+ϕTR-r

where λ, ϕ ∈ ℝ3.

The system we obtain after setting ∇Hc = 0 as Leeghim (2013) is:

fc=ηc = 0,R  - r= 0,ϕT  RX-rx+rµRt-rt= 0,J Δv1+1rϕTRxη3 Δv1-ϕTr Δv1=0,ϕT1rmur1η3x1+η3t1rx-mur1rx1+rt1=0. (6)

Fig. 3 Particle Swarm Optimization. 

5. SOLUTION METHOD

To find the roots of the systems (5) and (6) we use Newton’s multivariate method (Bate 1971; Leeghim 2013). That is, we want to calculate f (x) = 0 using

yn=yn-1-Jfyn-1-1fyn-1,

where Jfis the Jacobian matrix of f. The variables of the first problem are ys = (X, x, t, Δv0 {ϕ}) ∈ ℝ9 and the second one has yc = (X, x1, x, t1, t, Δv1 Δϕ) ∈ ℝ11.

For each case presented in § 6, it was necessary to find a seed (an initial value for Newton’s method). However, this seed must be very close to the solution for the method to converge (Local Convergence Theorem). Given the importance of calculating the seed, we used the following heuristic method to ap-proximate the solutions, which gave us very good results.

5.1. Particle Swarm Optimization

Particle Swarm Optimization (PSO) can be found in Kennedy et al. (1995); Clerc (2002); Parsopoulos & Vrahatis (2002); Conway (2010); Hvass (2010); Geetha et al. (2013). Let J : D ⊆ Rm → R be the function to be optimized, then:

  • N particles are randomly selected xi = (xi1,…,xim) ∈ D and also their initial velocities vi = (vi1,…,vim ) ∈ [0 1]m, i = 1,…,N.

Then, for each iteration k:

1. We find the minimum of {J(xj i)}j≤k (the minimum in the history of the specific particle) and it is set as pk i= xk i minfor each i = 1 N . Then we look for the minimum of {J(xm i)}m≤k i = 1 N (the minimum in the history of the entire set of particles) and after gk =xm imin.

2. The velocity vector is updated:

vik+1k+1=c1r1vvik+c2r2pik-xik+c3r3gik-gik,

where c1, c2, c3, ∈ (0 2] are adjustment parameters and r1 r 2 r 3 ∈ [0 1] are random numbers.

3. The particles are moved to their new position:

xik+1=xxik+vik+1.

To solve (3) or (4) we minimized Je = ||f || with the PSO, looking for different points xpso in ℝ9 or ℝ11. If Je ≈ 0, then we have a small region where the solution of f(x) = 0 is expected to be found. The region is defined by constructing an hyper-cube around the point found by the PSO, and then a grid of the region is made until Newton’s method converges to xnw for any of those seeds. Thus the desired solution is obtained.

Table 1 shows the efficiency of the PSO in the search for critical points in large regions (such as R11). Simulations were made using C and figures were obtained using Inkscape and MatLab, with a computer Lenovo ThinkCentre E73z i5-4430s, 2.7 GHz, with Windows 7.

TABLE 1 RESULTS OF PSO 

Processing Methods
Case Computing Time (s) J e(xpso ) # Iter (pso) Je( xnw) # Iter (nw) ||xpsoxnw||
LEO 25,94 10−4 5538 10−12 3641 1 61e−5
Mercury 2,62 0,3 241 10−12 60 2,28
Venus 1,05 0,1 313 10−12 95 0,38
Mars 0,67 0,5 72 10−12 61 0,47
Jupiter 2,63 0,5 321 10−12 36 0,71
Saturn 12,92 0,5 1804 10−12 26 0,27

6. RESULTS

The methods were applied to the following problem, proposed in Vallado (1997, p. 352), which deals with orbit transfers near Earth (Low Earth Orbit, LEO).

The Hubble telescope will be released from the space shuttle, which is in a circular orbit at 590 km from the Earth’s surface. The relative ejection speed (viewed from the shuttle) is [−0 1 −0 4 −0 2]T m/s. After 4 minutes, the Hubble needs to meet the shuttle. The change of RSW to IJK coordinates is shown in Vallado (1997, p. 367).

After applying the algorithm, the optimal orbit transfer without wait time has the parameters of the first row of Table 2 (for LEO), and has universal variables X = 136 206 and x = 136 207. With a wait time (shown by an asterisk in Table 2), X = 651 694, x1 = 265 476 and x = 386 215 are obtained. The orbits are shown in Figure 4.

TABLE 2 APPLICATIONS 

Change of Velocity
LEO t (min) | Δ v| (m/s) Δ vî Δ vĵ Δ vk^ Ignition
Leeghim 25.06 0.34 0.2 −0 27 −0 05 0
Leeghim* 71.04 0.15 0.0023 0 15 −0 01 48 min
Hohmann 48.24 0.2 · · · · · · · · · 48 min
Mercury t (d) | Δ v| (km/s) Δ vî Δ vĵ Δ vk^ Ignition
Leeghim 102 6.99 6.49 2.06 −1 61 01-01-2014
Leeghim* 99 6.96 6.53 1.8 −1 61 04-01-2014
Kemble 158 9.37 · · · · · · · · · 11-05-2012
Venus t (d) | Δ v| (km/s) Δ vî Δ vĵ Δ vk^ Ignition
Leeghim 85 20.33 20 24 −1 87 −0 05 01-01-2014
Leeghim* 142 2,42 −2 41 0 1 0 05 08-06-2014
Kemble 158 2.77 · · · · · · · · · 02-11-2013
Mars t (d) | Δ v| (km/s) Δ vî Δ vĵ Δ vk^ Ignition
Leeghim 200 3.15 −2 63 −1 13 1 25 01-01-2014
Leeghim* 207 2.99 −2 68 −0 63 −1 15 10-01-2014
Kemble 207 3.82 · · · · · · · · · 18-01-2014
Jupiter t (y) | Δ v| (km/s) Δ vî Δ vĵ Δ vk^ Ignition
Leeghim 7.19 9.69 −9 46 −1 28 0.99 01-01-2014
Leeghim* 2.22 9.39 −6 67 6 05 2.68 15-11-2016
Kemble 2.13 9.23 · · · · · · · · · 30-04-2009
Saturn t (y) | Δ v| (km/s) Δ vî Δ vĵ Δ vk^ Ignition
Leeghim 2.73 70.4 70 07 −6 34 −2 13 01-01-2014
Leeghim* 6.21 10.25 −9 33 −4 23 −0 02 13-01-2014
Kemble 9.2 10.49 · · · · · · · · · 22-12-2009

Fig 4. Transfers between Hubble (blue) - Shuttle (black). The color figure can be viewed online. 

The Hohmann transfer gives the optimal change of velocity in planar and circular orbits. Leeghim’s method deals with more general orbits: elliptic and hyperbolic; this method gives a smaller flight time.

It is also seen in Table 2 that for LEO orbits the time of flight of the transfer with wait time decreases 32 1% with respect to the one of Hohmann. In addition, Hohmann’s change in velocity is 33 3% higher than that with wait time, the latter evidencing significant fuel savings.

The following examples are the orbit transfers from Earth to the other planets of the solar system with a launch date on January 1st, 2014.

Figure 5 shows the current positions of Earth and Mercury for the assumed launch date, and the subsequent flight times. The resulting wait time is 3 days; with this, the change of velocity is reduced by one third and the flight time is reduced by 3 days with respect to the calculation without wait time.

Fig. 5 Transfers between Earth (blue) - Mercury (black). The color figure can be viewed online. 

These results are better than the values obtained in Kemble (2006, p. 50), where different launch dates for the optimization were assumed.

Computations were done without taking into account technical constraints. For example, the fastest space probe that NASA has launched is New Horizons with a |Δv| of 16 26 km/s relative to Earth, so the first |Δv| obtained without wait time is not feasible for such a mission, even though the flight time is the shortest.

Sometimes the direction of the Earth’s velocity [−29 61 − 6 41 0]T for January 1st, 2014) and the initial positions of the planets are not suitable for orbit transfers. This means that, in some cases, a very large change of velocity (which translates to a large fuel consumption) is required, as presented in the first row of Table 2 (for Venus). See Figure 6.

Fig. 6 Transfers between Earth (blue) - Venus (black). The color figure can be viewed online. 

Adding the wait time (167 days) results in a much lower change of velocity (12% of the first one, which requires less fuel), and the flight time increases by 57%, improving the results presented in Kemble (2006, p. 51) for comparable parameters.

The transfer of Earth-Mars orbit was calculated with wait time and without wait time. The difference lies in the fact that, with a wait time, the change of velocity was better and the flight time increased by a week. The results obtained are better than those of Kemble (2006, p. 53). See Figure 7.

Fig. 7 Transfers between Earth (blue) - Mars (black). The color figure can be viewed online. 

The system proposed in § 5.1 can have different solutions (local maximums of the index of performance J). The algorithm used to solve it finds many of them; here the optimal are shown.

For the example Earth-Jupiter, Figure 8, the three changes in velocity without a wait time, with wait time and as found in Kemble (2006, p. 55) are very similar, but the time of flight with wait time is smaller than the one without, and it is very close to the values obtained in Kemble (2006, p. 55). In this case, the eccentricity of the orbit transfers found is greater than the one of the previous cases (e = 0 75).

Fig. 8 Transfers between Earth (blue) - Jupiter (black). The color figure can be viewed online. 

The solution without a wait time for the system Earth-Saturn, Figure 9, is on a parabola (with eccentricity equal to 1, an uncommon result) and |Δ v| is extremely large. The transfer with wait time has an eccentricity of 0.82, giving a better change of velocity than the one presented in Kemble (2006, p. 56) and the flight time is 3 years shorter.

Fig. 9 Transfers between Earth (blue) - Saturn (black). The color figure can be viewed online. 

7. CONCLUSIONS

The two techniques presented here make it possible to optimize the fuel consumption. However, its use is purely theoretical, leaving aside the technical constraints, but considering only constraints of the trajectory to make a more general model. For Newton’s multivariate method, the search of the starting point was performed with the help of the PSO with excellent results. It perfectly bounds the search region of the starting point. The PSO is an heuristic method that needs no derivatives and has a rapid convergence. As future work we propose to use this method it in other situations and to perform the con-vergence analysis (local and semi-local) of the methods.

The authors are grateful with the young researcher fellowship of the Sergio Arboleda University from which this paper is one of its results.

REFERENCES

Bate, R. R., Mueller, D. D., & White, J. E. 1971, Fundamentals of Astrodynamics (New York: Dover Publications) [ Links ]

Clerc, M. & Kennedy, J. 2002, ITEVC, 6(1), 58 [ Links ]

Conway, B. 2010, Spacecraft Trajectory Optimization (Cambridge, MA: CUP) [ Links ]

Danby, J. M. A. 1992, Fundamentals of Celestial Mechanics, (New York, NY: Macmillan) [ Links ]

Elipe, A., Montijano, J. I., Rández, L., & Calvo, M. 2017, CeMDA, 129, 415 [ Links ]

Gang, Z., Xiangyu, Z., & Xibin, C. 2014, ChJA, 27(3), 577 [ Links ]

Geetha, T., Sowmiya, B., JayaKumar, L., & Sarang Sukumar, A. 2013, IJETAE, 3(4), 672 [ Links ]

Hvass, M. 2010, Tuning And Simplifying Heuristical Optimization [Thesis], (Southampton: University of Southampton) [ Links ]

Kemble, S. 2006, Interplanetary Missions Analysis and Design, (Berlin: Springer-Verlag Berlin Heidelberg) [ Links ]

Kennedy, J. & Eberhart, R. 1995, Particle Swarm Optimization, (Perth, WA: IEEEP) [ Links ]

Leeghim, H. 2013, CeMDA , 115, 1 [ Links ]

Parsopoulos, K. & Vrahatis, M. 2002, In Frontiers in Artificial Intelligence and Applications, ed. P. Sinkad et al. (Amsterdam: IOS Press), 76, 214 [ Links ]

Sharaf, M. A. & Saad, A. S. 2016, RMxAA, 52, 283 [ Links ]

Shan, J. & Ren, Y. 2014, Aerospace Science and Technology, 36, 114 [ Links ]

Vallado, D. A. 1997, Fundamentals of Astrodynamics and Applications (New York, NY: McGraw-Hill) [ Links ]

Villanueva, Y. 2013, Lagrangian Optimization in Inter-planetary Trajectories [Thesis], (Bogotá: Universidad Sergio Arboleda) [ Links ]

Zotes, F. & Santos, M. 2012, Engineering Applications of Artificial Intelligence, 25(1), 189 [ Links ]

Received: April 24, 2019; Accepted: May 09, 2019

L. M. Echeverry and Y. Villanueva: School of Exact Sciences and Engineering, Sergio Arboleda University, Bogotá, Colombia (lechever@uniandes.edu.co, yovaniing@gmail.com).

L. M. Echeverry: Department of Mathematics, University of the Andes, Bogotá, Colombia (lechever@uniandes.edu.co).

Y. Villanueva: Institute of Mathematics and Statistics, Federal University of Goiás, Goiania, Brazil (yovaniing@gmail.com).

Creative Commons License Este es un artículo publicado en acceso abierto bajo una licencia Creative Commons