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.
Loading packages¶
First we load some python packages for plotting and linear algebra.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
Then we load godot specific libraries and the pygmo optimisation package
from godot.core import tempo, util, constants
from godot.model import common
from godot import cosmos
import pygmo
Universe¶
Load the universe configuration and create the universe object. This is a barebones configuration with using JPL ephemeris for planets and some planetary axes models. The dynamics include Sun, Earth and Mars, and the spacecraft SEP engine is configured with fixed values for Isp and thrust.
util.suppressLogger()
uni_config = cosmos.util.load_yaml("universe.yml")
uni = cosmos.Universe(uni_config)
Trajectory¶
Load the trajectory configuration and create the trajectory object using the universe object.
The trajectory optimisation problem we are solving is an Earth to Mars transfer with solar electric propulsion (SEP). It contains two segments (separate propagations), one forward from Earth and one backward from Mars. There is a matching point in between. Each arc contains a continuous thrust arc with mandatory coast after launch and before arrival (30 days).
tra_config = cosmos.util.load_yaml("trajectory.yml")
tra = cosmos.Trajectory(uni, tra_config)
Plotting¶
We generate a plotting function to create visualization for our trajectory in the solar system.
def man_data(name):
man = tra.getManoeuvreBook()[name]
start = tra.point(man.start)
end = tra.point(man.end)
ran = tempo.EpochRange(start, end)
grid = ran.createGrid(1.0 * tempo.SecondsInDay)
mass = uni.evaluables.get("SC_mass")
t = [e.mjd() for e in grid]
r = np.asarray([uni.frames.vector3("Sun", "SC_center", "EMC", e) / constants.AU for e in grid])
man.model.setActive(True)
thr = np.asarray([np.linalg.norm(man.model.eval(e)) * mass.eval(e) for e in grid])
man.model.setActive(False)
return t, r, thr
def plot():
tra.compute(False)
ran = tempo.EpochRange(tempo.Epoch("2020-01-01 TDB"), tempo.Epoch("2022-01-01 TDB"))
earth = np.asarray([uni.frames.vector3("Sun", "Earth", "EMC", e) / constants.AU for e in ran.createGrid(1.0 * tempo.SecondsInDay)])
mars = np.asarray([uni.frames.vector3("Sun", "Mars", "EMC", e) / constants.AU for e in ran.createGrid(1.0 * tempo.SecondsInDay)])
def pos(e):
try:
return uni.frames.vector3("Sun", "SC_center", "EMC", e) / constants.AU
except:
return [np.nan] * 3
grid = tra.range().createGrid(1.0 * tempo.SecondsInDay)
t = [e.mjd() for e in grid]
thr = np.asarray([uni.evaluables.get("SC_sep_thrust").eval(e) for e in grid])
x = np.asarray([pos(e) for e in grid])
man1_t, man1_r, man1_thr = man_data("man1")
man2_t, man2_r, man2_thr = man_data("man2")
match1_ep = tra.point("match_inter_left")
match1_val = uni.evaluables.get("SC_sep_thrust").eval(match1_ep)
match2_ep = tra.point("match_inter_right")
match2_val = uni.evaluables.get("SC_sep_thrust").eval(match2_ep)
plt.figure(figsize=(8, 8))
plt.xlabel("EMC X (AU)")
plt.ylabel("EMC Y (AU)")
plt.plot(earth[:, 0], earth[:, 1], "--", label="Earth")
plt.plot(mars[:, 0], mars[:, 1], "--", label="Mars")
plt.plot(x[:, 0], x[:, 1], "-k", linewidth=2, label="Transfer")
plt.plot(man1_r[:, 0], man1_r[:, 1], "-r", linewidth=4, label="SEP 1")
plt.plot(man2_r[:, 0], man2_r[:, 1], "-g", linewidth=4, label="SEP 2")
plt.plot(x[0, 0], x[0, 1], "ok", ms=10)
plt.plot(x[-1, 0], x[-1, 1], "ok", ms=10)
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.legend()
plt.grid()
plt.figure(figsize=(8, 4))
plt.xlabel("Time (MJD)")
plt.ylabel("SEP thrust (mN)")
plt.plot(t, 1e6 * thr, "--k", linewidth=1, label="Max thrust")
plt.plot([tra.range().start().mjd()], 1e6 * thr[0], "o", ms=10, label="Earth")
plt.plot([tra.range().end().mjd()], 1e6 * thr[-1], "o", ms=10, label="Mars")
plt.plot([match1_ep.mjd()], 1e6 * match1_val, "dk", ms=5, label="Match")
plt.plot([match2_ep.mjd()], 1e6 * match2_val, "dk", ms=5)
plt.plot(man1_t, 1e6 * man1_thr, "-r", linewidth=3, label="SEP 1")
plt.plot(man2_t, 1e6 * man2_thr, "-g", linewidth=3, label="SEP 2")
plt.xlim([tra.range().start().mjd(), tra.range().end().mjd()])
plt.ylim([0, 1000])
plt.legend()
plt.grid()
We compute the initial guess of the transfer trajectory and visualize it. Immediately we can see that there is a large discontinuity at the matching point and this optimisation problem needs to be solved.
plot()
Summary¶
We add a function to summarize some results.
def summary():
print(f"Escape velocity: {uni.parameters.get('Departure_SC_center_vin').getPhysicalValue():.3f} km/s")
print(f"Launch mass: {uni.parameters.get('Departure_SC_mass').getPhysicalValue():.1f} kg")
print(f"Arrival mass: {uni.parameters.get('Arrival_SC_mass').getPhysicalValue():.1f} kg")
print(f"Total delta-V: {uni.parameters.get('Arrival_SC_dv').getPhysicalValue():.3f} km/s")
summary()
Problem¶
Load the problem configuration and create the problem object using the universe, a list of considered trajectories. The problem defines free variables and optionally their bounds/scales as well as additional constraints and the objective.
In this optimisation example we fix the departure/arrival dates and escape/arrival infinite velocities and only optimise the V-infinity directions, as well as the low-thrust maoneuvre directions.
pro_config = cosmos.util.load_yaml("problem.yml")
pro = cosmos.Problem(uni, [tra], pro_config)
Now it is time to create the pygmo problem class that simply takes ours, and define some constraint tolerances
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
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
algo = pygmo.algorithm(pygmo.nlopt("slsqp"))
algo.set_verbosity(1)
pop = algo.evolve(pop)
We can visualize the final converged trajectory.
plot()
summary()
Finally, we can save the optimised and updated trajectory configuration
up = tra.applyParameterChanges()
tra_config = cosmos.util.deep_update(tra_config, up)
cosmos.util.save_yaml(tra_config, "trajectory_out.yml")
Custom thrust model¶
We will create a new universe to show how easy it is to implement a new thrust model.
uni_config = cosmos.util.load_yaml("universe.yml")
uni = cosmos.Universe(uni_config)
Our "complex" thrust model will replace the fixed thrust model with one that depends on the distance from the Sun. We will use 500 mN thrust available at 1 AU, and scale it with the inverse square of the distance from the Sun.
class ThrustModel(common.ScalarTimeEvaluable):
def __init__(self, uni, thrust_at_1AU):
super().__init__()
self.uni = uni
self.thrust_at_1AU = thrust_at_1AU
def eval(self, e):
dist = self.uni.frames.distance("Sun", "SC_center", e) / constants.AU
thr = self.thrust_at_1AU / (dist * dist)
return thr
thrust_at_1AU = 500e-6 # in mN
thrust_model = ThrustModel(uni, thrust_at_1AU)
uni.evaluables.replace("SC_sep_thrust", thrust_model)
As initial guess, we us the converged trajectory using the fixed thrust model.
tra_config = cosmos.util.load_yaml("trajectory_out.yml")
tra = cosmos.Trajectory(uni, tra_config)
We visualize the initial guess which uses the new thrust model. We expect to see some discontinuities as we updated the model.
plot()
Finally, we optimise the trajectory.
pro_config = cosmos.util.load_yaml("problem.yml")
pro = cosmos.Problem(uni, [tra], pro_config)
problem = pygmo.problem(pro)
conTol = 1e-6
problem.c_tol = [conTol] * problem.get_nc()
x0 = pro.get_x()
pop = pygmo.population(problem, 0)
pop.push_back(x0)
algo = pygmo.algorithm(pygmo.nlopt("slsqp"))
algo.set_verbosity(1)
pop = algo.evolve(pop)
And visualize the new solution with the inverse square thrust model.
plot()
summary()
We save this trajectory as well.
up = tra.applyParameterChanges()
tra_config = cosmos.util.deep_update(tra_config, up)
cosmos.util.save_yaml(tra_config, 'trajectory_out2.yml')
Using launcher performance model¶
In the following we further complicate the problem by replacing the fixed departure boundary conditions (escape velocity, mass) with one that links them through the launcher performance model.
We use a simple linear model in C3 as in most cases that is quite close to reality.
def launcher_model(vinf):
c3 = vinf * vinf
coef0 = 4000.0
coef1 = -250.0
mass = coef0 + coef1 * c3
return mass
vin = np.linspace(0, 3)
mass = [launcher_model(x) for x in vin]
plt.plot(vin, mass, lw=2)
plt.grid()
plt.xlabel("Escape velocity (km/s)")
plt.ylabel("Injected mass (kg)")
plt.show()
We recreate the same universe as before including the custom thrust model
uni_config = cosmos.util.load_yaml("universe.yml")
uni = cosmos.Universe(uni_config)
thrust_model = ThrustModel(uni, thrust_at_1AU)
uni.evaluables.replace("SC_sep_thrust", thrust_model)
We load the converged trajectory from the previous step.
tra_config = cosmos.util.load_yaml("trajectory_out2.yml")
tra = cosmos.Trajectory(uni, tra_config)
In addition, now we add an evaluable that defines the constraint defect of the launcher performance. We add it to the universe evaluables so that it is avaiable to be used later in the problem definition.
class LauncherConstraint(common.ScalarTimeEvaluable):
def __init__(self, uni):
super().__init__()
self.uni = uni
def eval(self, e):
mass = self.uni.parameters.get("Departure_SC_mass").eval(e)
vinf = self.uni.parameters.get("Departure_SC_center_vin").eval(e)
out = launcher_model(vinf) - mass
return out
launcher_con = LauncherConstraint(uni)
uni.evaluables.add("launcher_con", launcher_con)
And now we load a different problem configuration, that allows the initial escape velocity and mass to be optimised, as well as defining the additional constraint. Moreover, we need to switch the objective from minimising the delta-V to optimising the final mass, as now that is our final merit value.
pro_config = cosmos.util.load_yaml("problem2.yml")
pro = cosmos.Problem(uni, [tra], pro_config)
# We could also add constraint directly via the problem interface
# pro.addConstraint(1e-3 * launcher_con)
Carry out the new optimisation.
def optimise():
problem = pygmo.problem(pro)
conTol = 1e-6
problem.c_tol = [conTol] * problem.get_nc()
x0 = pro.get_x()
pop = pygmo.population(problem, 0)
pop.push_back(x0)
algo = pygmo.algorithm(pygmo.nlopt("slsqp"))
algo.set_verbosity(1)
pop = algo.evolve(pop)
optimise()
And plot the converged trajectory. We can see that the thrust arcs are not saturating the transfer, as the mass is pushed to its limit to optimise the final mass. This results in lower escape velocity and lower thrust-to-mass ratio, so the optimiser increases the thrust act durations to their limit.
plot()
summary()
Save resulting trajectory.
up = tra.applyParameterChanges()
tra_config = cosmos.util.deep_update(tra_config, up)
cosmos.util.save_yaml(tra_config, "trajectory_out3.yml")
Mass dependent thrust¶
In the final iteration of our example, we will include a system design decision. The thrust available at 1 AU will depend on the initial wet mass of the SC. This can be thought as optimising the number of solar arrays but taking into account their additional mass. In most cases this decision variable would be discrete, but here we just use a continuous variable for simplicity.
We define the model that computes the avaialble thrust at 1 AU from the wet mass. We include some 2nd order term to account for additional dry mass needed to support the more bulky power infrastructure.
def thrust_at_1AU_from_mass(mass):
m0 = 2000.0
dm = mass - m0
coef1 = 0.5e-6
coef2 = -8e-11
return coef1 * dm + coef2 * dm * dm
mass = np.linspace(2000, 5000)
thr = np.asarray([thrust_at_1AU_from_mass(x) for x in mass])
plt.plot(mass, 1e6 * thr, lw=2)
plt.grid()
plt.xlabel("Wet mass (kg)")
plt.ylabel("Thrust (mN)")
plt.plot()
plt.show()
No we create the thrust model that uses this new assumption.
class ComplexThrustModel(common.ScalarTimeEvaluable):
def __init__(self, uni):
super().__init__()
self.uni = uni
def eval(self, e):
dist = self.uni.frames.distance("Sun", "SC_center", e) / constants.AU
wet_mass = self.uni.parameters.get("Departure_SC_mass").eval(e)
thr_at_1AU = thrust_at_1AU_from_mass(wet_mass)
thr = thr_at_1AU / (dist * dist)
return thr
Similarly set up the universe as before.
uni_config = cosmos.util.load_yaml("universe.yml")
uni = cosmos.Universe(uni_config)
thrust_model = ComplexThrustModel(uni)
uni.evaluables.replace("SC_sep_thrust", thrust_model)
launcher_con = LauncherConstraint(uni)
uni.evaluables.add("launcher_con", launcher_con)
As well as loading the previous converged trajectory.
tra_config = cosmos.util.load_yaml("trajectory_out3.yml")
tra = cosmos.Trajectory(uni, tra_config)
And plotting the initial guess, noticing the discontinuity due to the different thrust model.
plot()
We can reuse the same problem configuration from before.
pro_config = cosmos.util.load_yaml('problem2.yml')
pro = cosmos.Problem(uni, [tra], pro_config)
And optimise the problem.
uni.parameters.get("match_inter_left_dt").setScaledValue(1e-9)
optimise()
Finally, plot the solution.
plot()
summary()
And save the final trajectory.
up = tra.applyParameterChanges()
tra_config = cosmos.util.deep_update(tra_config, up)
cosmos.util.save_yaml(tra_config, 'trajectory_out4.yml')