Trajectory optimisation example

This example shows how to optimise a trajectory using the GODOT cosmos module and the third-party open source optimisation package, pagmo/pygmo.

Pagmo/Pygmo is developed by the Advanced Concept Team (ACT) in ESTEC, and contains many freely available global and local optimisation algorithms.

First we load some python packages for plotting and linear algebra.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Then we load godot specific libraries and the pygmo optimisation package

In [2]:
from godot.core import tempo
from godot import cosmos
import pygmo

# optionally avoid verbose logging messages
import godot.core.util as util
util.suppressLogger()

Load the universe configuration and create the universe object

In [3]:
uniConfig = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uniConfig)

Load the trajectory configuration and create the trajectory object using the universe object

In [4]:
traConfig = cosmos.util.load_yaml('trajectory.yml')
tra = cosmos.Trajectory(uni, traConfig)

Load the problem configuration and create the problem object using the universe, a list of considered trajectories

In [5]:
proConfig = cosmos.util.load_yaml('problem.yml')
pro = cosmos.Problem(uni, [tra], proConfig)

Now it is time to create the pygmo problem class that simply takes ours, and define some constraint tolerances

In [6]:
problem = pygmo.problem(pro)
conTol = 1e-6
problem.c_tol = [conTol] * problem.get_nc()

We create the pygmo population to be optimised by grabbing the initial value of the free parameters

In [7]:
x0 = pro.get_x()
pop = pygmo.population(problem, 0)
pop.push_back(x0)

Finally we use the NLOPT open source package to optimise the problem with the well-known SLSQP algorithm

In [8]:
algo = pygmo.algorithm(pygmo.nlopt("slsqp"))
algo.set_verbosity(1)
pop = algo.evolve(pop)

After evaluation the trajectory, we can request the solutions of specific timeline elements

In [9]:
sol = tra.getTimelineSolution()
for entry in sol:
    for item in entry:
        print('%24s%15.6f' % (item.name, item.epoch.mjd()))

                  launch       0.000000
              man1_start       0.426221
                man1_end       0.532101
             match1_left       1.000000
            match1_right       1.000000
                man2_end       1.533049
              man2_start       1.568330
                 arrival       2.000000

For each propagation arc, we generate some plotting data for the position vector, propagated mass and accumulated delta-V

In [10]:
data = list()
for entry in sol:
    for left, right in zip(entry, entry[1:]):
        grid = tempo.EpochRange(left.epoch, right.epoch).createGrid(1e-3 * 86400)
        t = [e.mjd() for e in grid]
        x = np.asarray([uni.frames.vector3('Earth', 'GeoSat_center', 'ICRF', e) for e in grid])
        m = [uni.evaluables.get('GeoSat_mass').eval(e) for e in grid]
        v = [uni.evaluables.get('GeoSat_dv').eval(e) for e in grid]
        data.append((t, x, m, v))

Let's plot the X-Y components of the evolution of the position vector

In [11]:
plt.grid()
plt.xlabel('X pojection [km]')
plt.ylabel('Y pojection [km]')
for i, entry in enumerate(data):
    c = [u'b', u'g', u'r', u'c', u'm', u'y', u'k'][i]
    plt.plot([entry[1][0, 0]], [entry[1][0, 1]], 'o', color=c)
    plt.plot(entry[1][:, 0], entry[1][:, 1], '-', color=c)
    plt.plot([entry[1][-1, 0]], [entry[1][-1, 1]], 'o', color=c)

Let's plot the 3D Earth-centered trajectory

In [12]:
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.grid()
for i, entry in enumerate(data):
    c = [u'b', u'g', u'r', u'c', u'm', u'y', u'k'][i]
    ax.plot([entry[1][0, 0]], [entry[1][0, 1]], [entry[1][0, 2]], 'o', color=c)
    ax.plot(entry[1][:, 0], entry[1][:, 1], entry[1][:, 2], '-', color=c)
    ax.plot([entry[1][-1, 0]], [entry[1][-1, 1]], [entry[1][-1, 2]], 'o', color=c)

Let's plot the mass evolution

In [13]:
plt.grid()
plt.xlabel('Time [mjd]')
plt.ylabel('Mass [kg]')
for i, entry in enumerate(data):
    c = [u'b', u'g', u'r', u'c', u'm', u'y', u'k'][i]
    plt.plot(entry[0][0], entry[2][0], 'o', color=c)
    plt.plot(entry[0], entry[2], '-', color=c)
    plt.plot(entry[0][-1], entry[2][-1], 'o', color=c)

Let's plot the accumulated delta-V evolution

In [14]:
plt.grid()
plt.xlabel('Time [mjd]')
plt.ylabel('Delta-V [km/s]')
for i, entry in enumerate(data):
    c = [u'b', u'g', u'r', u'c', u'm', u'y', u'k'][i]
    plt.plot(entry[0][0], entry[3][0], 'o', color=c)
    plt.plot(entry[0], entry[3], '-', color=c)
    plt.plot(entry[0][-1], entry[3][-1], 'o', color=c)

Finally, we can save the optimised and updated trajectory configuration

In [15]:
up = tra.applyParameterChanges()
traConfig = cosmos.util.deep_update(traConfig, up)
cosmos.util.save_yaml(traConfig, 'trajectory_out.tra.yml')
In [ ]: