Multiple Gravity Assist (MGA)

In this tutorial, we use the pykep.trajopt.mga to design an interplanetary trajectory. The MGA encoding of an interplanetary trajectory models multiple (pre-defined) planetary encounters, and enables for propulsion manovures in the form of instantaneous \(\Delta V\) immediately after each fly-by. The encoding is useful as a preliminar mission design tool automating the search for good planetary geometries giving the possibility to design trajectories through a process mimicking natural evolution.

This encoding has also been used in the past to define instances of real-world optimization problems to be used in benchmarking optimization tools suitable for interplanetary mission design.

We start, as often, with some fundamental imports:

import pykep as pk
import pygmo as pg
import matplotlib.pyplot as plt

We will study the possibility of establishing an Earth return trajectory visiting venus in two different epochs and mars in between:

earth = pk.planet(pk.udpla.jpl_lp("earth"))
venus = pk.planet(pk.udpla.jpl_lp("venus"))
mars = pk.planet(pk.udpla.jpl_lp("mars"))
udp = pk.trajopt.mga(
    seq=[
        earth,
        venus,
        mars,
        venus,
        earth
    ],
    tof_encoding = "direct",
    t0=[0, 1000],
    tof=[[30, 200], [30, 300], [30, 350], [30, 300]],
    vinf=5.,
)

To solve this problem we will use an adaptive version of differential evolution, even if other meta-heuristic would also be possible. We run the algorithm in mult-start, to get a sense of the fitness landscape and complexity.

prob = pg.problem(udp)
# CMA-ES is also possible
#uda = pg.cmaes(1500, force_bounds=True, sigma0=0.5, ftol=1e-4)
# But we prefer a self adaptive version of differential evolution here
uda = pg.sade(2500, ftol=1e-4, xtol=1e-4)
algo = pg.algorithm(uda)
res = list()
for i in range(10):
    pop = pg.population(prob, 20)
    pop = algo.evolve(pop)
    res.append([pop.champion_f, pop.champion_x])
    print(i, pop.champion_f[0], flush=True)
    
best_x = sorted(res, key =  lambda x: x[0][0])[0][1]

0 3783.1332872905305

1 3787.606982393975

2 3784.329725783436

3 3785.0187280829564

4 3784.2894750411438

5 3780.655789075206

6 3786.283132277777

7 3784.694024584519

8 3787.39782021129

9 3783.83711704293

Plotting the trajectory

There are several ways one can plot the solution, depending on what we want to ‘see’. Lets start with one “at-a-glance” plot:

ax = udp.plot(best_x, figsize=(7, 7))
ax.view_init(90, 0)
ax.axis("off");

If we need to disentangle the sequence of the varioous legs we may also subdivide them in subplots:

# The figure
fig = plt.figure(figsize=(18,18))

# The single leg plots
for i in range(1,5):
    ax = fig.add_subplot(1, 5, i, projection="3d")
    ax.view_init(90, 0)
    ax = udp.plot(best_x, ax=ax, leg_ids=[i-1])
    ax.axis("off");
    
# All the legs
ax = fig.add_subplot(1, 5, 5, projection="3d")
ax.view_init(90, 0)
ax = udp.plot(pop.champion_x, ax=ax)

… or show the various xy, xy, xz views as well as a full 3D one:

fig = plt.figure(figsize=(8,8))

ax = fig.add_subplot(2, 2, 1, projection="3d")
ax = udp.plot(best_x, ax=ax)
ax.view_init(45, 30)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
ax.get_zaxis().set_ticks([])

ax = fig.add_subplot(2, 2, 2, projection="3d")
ax = udp.plot(best_x, ax=ax)
ax.view_init(90, 0)
ax.axis("off")

ax = fig.add_subplot(2, 2, 3, projection="3d")
ax = udp.plot(best_x, ax=ax)
ax.view_init(0, 0)
ax.axis("off")

ax = fig.add_subplot(2, 2, 4, projection="3d")
ax = udp.plot(best_x, ax=ax)
ax.view_init(0, 90)
ax.axis("off");

Analyzing the trajectory

For a quick analysis of the trajectory we may use the pretty method.

ax = udp.pretty(best_x)

Multiple Gravity Assist (MGA) problem: 
	Planet sequence:  ['earth - jpl_lp', 'venus - jpl_lp', 'mars - jpl_lp', 'venus - jpl_lp', 'earth - jpl_lp']
	Encoding for tofs:  direct
	Orbit Insertion:  False
Departure:  earth - jpl_lp
	Epoch:  933.511649745537  [mjd2000]
	Spacecraft velocity:  [23547.740944388835, 12502.609390978907, 3065.7471491197957] [m/s]
	Hyperbolic velocity:  4147.943728828037 [m/s]
	Initial DV:  0 [m/s]
Fly-by:  venus - jpl_lp
	Epoch:  1076.3371094871202  [mjd2000]
	DV:  0.009356465452583507 [m/s]
Fly-by:  mars - jpl_lp
	Epoch:  1233.7768243131907  [mjd2000]
	DV:  1.5537134459009394e-06 [m/s]
Fly-by:  venus - jpl_lp
	Epoch:  1534.5114982238176  [mjd2000]
	DV:  0.00039539364297525026 [m/s]
Arrival:  earth - jpl_lp
	Epoch:  1691.0412162453788  [mjd2000]
	Spacecraft velocity:  [15614.999213034618, 21652.321501762814, 2460.4862843520373] [m/s]
	Arrival DV:  3780.646035662397 [m/s]
Time of flights:  [142.82545974 157.43971483 300.73467391 156.52971802] [days]

Total DV:  3780.655789075206