Global optimization of an MGA low-thrust trajectory¶

The example demonstrates the use of pykep/pygmo to perform global optimization of a multiple leg interplanetary trajectory over large launch windows. In particular, it defines a transfer between Earth and Mercury making use of a Venus intermediate fy-by as an global optimization problem (using the module trajopt udps) and it then attempts to find one solution using the Monotonic Basin Hopping meta-algorithm connected to an SQP local optimization technique.

This solution technique (MBH + SQP) was developed by members of the pykep team and the original publication can be found here:

Yam, C. H., D. D. Lorenzo, and D. Izzo. “Low-thrust trajectory design as a constrained global optimization problem.” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering 225, no. 11 (2011): 1243-1251.

A., Cassioli, D. Izzo, Di Lorenzo, D., Locatelli, M. and Schoen, F.: “Global optimization approaches for optimal trajectory planning.” In Modeling and optimization in space engineering, pp. 111-140. Springer, New York, NY, 2012.

[1]:

# Imports
import pykep as pk
import pygmo as pg
import pygmo_plugins_nonfree as ppnf
import numpy as np
from pykep.examples import add_gradient

# Plotting imports
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-d4897bb3ff63> in <module>
      2 import pykep as pk
      3 import pygmo as pg
----> 4 import pygmo_plugins_nonfree as ppnf
      5 import numpy as np
      6 from pykep.examples import add_gradient

ModuleNotFoundError: No module named 'pygmo_plugins_nonfree'

[2]:

# We define the optimization problem udp (User Defined Problem): an Earth-Venus-Mercury transfer with low-thrust
udp = pk.trajopt.mga_lt_nep(
     seq = [pk.planet.jpl_lp('earth'), pk.planet.jpl_lp('venus'), pk.planet.jpl_lp('mercury')],
     n_seg = [5, 20],
     t0 = [3000, 4000], # This is in mjd2000
     tof = [[100, 1000], [200, 2000]], # This is in days
     vinf_dep = 3., # This is in km/s
     vinf_arr = 2., # This is in km/s
     mass = [1000., 2000.0],
     Tmax = 0.5,
     Isp = 3500.0,
     fb_rel_vel = 6., # This is in km/s
     multi_objective = False,
     high_fidelity = False
)
prob = pg.problem(udp)
prob.c_tol = 1e-4
print(prob)

Problem name: <class 'pykep.trajopt._mga_lt_nep.mga_lt_nep'>
        C++ class name: pybind11::object

        Global dimension:                       92
        Integer dimension:                      0
        Fitness dimension:                      44
        Number of objectives:                   1
        Equality constraints dimension:         15
        Inequality constraints dimension:       28
        Tolerances on constraints: [0.0001, 0.0001, 0.0001, 0.0001, 0.0001, ... ]
        Lower bounds: [3000, 100, 1000, -3000, -3000, ... ]
        Upper bounds: [4000, 1000, 2000, 3000, 3000, ... ]
        Has batch fitness evaluation: false

        Has gradient: false
        User implemented gradient sparsity: false
        Has hessians: false
        User implemented hessians sparsity: false

        Fitness evaluations: 0

        Thread safety: none

[8]:

# We define the optimization algorithm uda (User Defined Algorithm). In this case we use SNOPT from the
# module pgmo_plugins_nonfree. Note that we point to the snopt7_c library in our system.
# In such a library is not available, nlopt sqlsqp can also be used here ....
uda = ppnf.snopt7(screen_output = False, library = "/usr/local/lib/libsnopt72_c.so", minor_version = 2)
uda.set_integer_option("Major iterations limit", 1000)
uda.set_integer_option("Iterations limit", 200000)
uda.set_numeric_option("Major optimality tolerance", 1e-2)
uda.set_numeric_option("Major feasibility tolerance", 1e-8)
uda2 = pg.mbh(uda, 3, 0.05)
algo = pg.algorithm(uda2)
algo.set_verbosity(1)

[9]:

# We create a population of 100 random initial guesses
pop = pg.population(prob, 100)
# And optimize
pop = algo.evolve(pop)

[15]:

print("Is feasible: ", prob.feasibility_f(pop.champion_f))
print("Final Mass at Mercury: ", -pop.champion_f[0], "Kg")

Is feasible:  True
Final Mass at Mercury:  1227.9619205547747 Kg

[14]:

# We plot
mpl.rcParams['legend.fontsize'] = 10

# Create the figure and axis
fig = plt.figure(figsize = (16,5))
ax1 = fig.add_subplot(1, 3, 1, projection='3d')
udp.plot(pop.champion_x, axes = ax1)

ax2 = fig.add_subplot(1, 3, 2, projection='3d')
ax2.view_init(90, 0)
udp.plot(pop.champion_x, axes = ax2)

ax3 = fig.add_subplot(1, 3, 3, projection='3d')
ax3.view_init(0,0)
udp.plot(pop.champion_x, axes = ax3)

[14]:

<matplotlib.axes._subplots.Axes3DSubplot at 0x7f5327950470>

[ ]: