Planet to planet low-thrust#
In this tutorial we show the use of the pykep.trajopt.direct_pl2pl to find a low-thrust trajectory connecting two moving planets.
The decision vector for this class, compatible with pygmo [BI20] UDPs (User Defined Problems), is:
containing the starting epoch \(t_0\) as a MJD2000, the final mass \(m_f\) as well as the starting and final \(V^{\infty}\), throttles and the time-of-flight \(T_{tof}\).
Note
This notebook makes use of the commercial solver SNOPT 7 and to run needs a valid snopt_7_c library installed in the system. In case SNOPT7 is not available, you can still run the notebook using, for example uda = pg.algorithm.nlopt("slsqp") with minor modifications.
Basic imports:
import pykep as pk
import numpy as np
import time
import pygmo as pg
import pygmo_plugins_nonfree as ppnf
import time
from matplotlib import pyplot as plt
We start defining the problem data.
# Problem data
mu = pk.MU_SUN
max_thrust = 0.3
isp = 3000
# Source and destination planets
earth = pk.planet_to_keplerian(
pk.planet(pk.udpla.jpl_lp(body="EARTH")), when=pk.epoch(5000)
)
mars = pk.planet_to_keplerian(
pk.planet(pk.udpla.jpl_lp(body="MARS")), when=pk.epoch(5000)
)
# Initial state
ms = 1500.0
# Number of segments
nseg = 10
We here instantiate two different versions of the same UDP (User Defined Problem), with analytical gradients and without.
For the purpose of this simple notebook we choose a relatively simple Earth to Mars transfer with an initial \(V_{\infty}\) of 3 km/s.
udp_nog = pk.trajopt.direct_pl2pl(
pls=earth,
plf=mars,
ms=ms,
mu=pk.MU_SUN,
max_thrust=max_thrust,
isp=3000,
t0_bounds=[6700.0, 6800.0],
tof_bounds=[200.0, 300.0],
mf_bounds=[1300., ms],
vinfs=3.,
vinff=0.,
nseg=nseg,
cut=0.6,
mass_scaling=ms,
r_scaling=pk.AU,
v_scaling=pk.EARTH_VELOCITY,
with_gradient=False,
)
udp_g = pk.trajopt.direct_pl2pl(
pls=earth,
plf=mars,
ms=ms,
mu=pk.MU_SUN,
max_thrust=max_thrust,
isp=3000,
t0_bounds=[6700.0, 6800.0],
tof_bounds=[200.0, 300.0],
mf_bounds=[1300., ms],
vinfs=3.,
vinff=0.,
nseg=nseg,
cut=0.6,
mass_scaling=ms,
r_scaling=pk.AU,
v_scaling=pk.EARTH_VELOCITY,
with_gradient=True,
)
Analytical performances of the analytical gradient#
And we take a quick look at the performances of the analytical gradient with respect to the numerically computed one.
# We need to generste a random chromosomes compatible with the UDP where to test the gradient.
prob_g = pg.problem(udp_g)
pop_g = pg.population(prob_g, 1)
First the analytical gradient:
%%timeit
udp_g.gradient(pop_g.champion_x)
166 μs ± 239 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
Then a simple numerical gradient based on finite differences:
%%timeit
pg.estimate_gradient(udp_g.fitness, pop_g.champion_x)
1.03 ms ± 618 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
Then a higher order numerical gradient:
%%timeit
pg.estimate_gradient_h(udp_g.fitness, pop_g.champion_x)
3.12 ms ± 34.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
The analytical gradient its exact and faster, seems like a no brainer to use it.
In reality, the effects on the optimization technique used are not straightforward, making the option to use numerical gradients still interesting in some, albeit rare, cases.
Solving the low-thrust transfer#
We define (again) the optimization problem, and set a tolerance for pagmo to be able to judge the relative value of two individuals.
Note
This tolerance has a different role from the numerical tolerance set in the particular algorithm chosen to solve the problem and is only used by the pagmo machinery to decide outside the optimizer whether the new proposed indivdual is better than what was the previous champion.
prob_g = pg.problem(udp_g)
prob_g.c_tol = 1e-6
… and we define an optimization algorithm.
snopt72 = "/Users/dario.izzo/opt/libsnopt7_c.dylib"
uda = ppnf.snopt7(library=snopt72, minor_version=2, screen_output=False)
uda.set_integer_option("Major iterations limit", 2000)
uda.set_integer_option("Iterations limit", 20000)
uda.set_numeric_option("Major optimality tolerance", 1e-3)
uda.set_numeric_option("Major feasibility tolerance", 1e-11)
#uda = pg.nlopt("slsqp")
algo = pg.algorithm(uda)
We solve the problem from random initial guess ten times and only save the result if a feasible solution is found (as defined by the criterias above)
masses = []
xs = []
for i in range(10):
pop_g = pg.population(prob_g, 1)
pop_g = algo.evolve(pop_g)
if(prob_g.feasibility_f(pop_g.champion_f)):
print(".", end="")
masses.append(pop_g.champion_x[1])
xs.append(pop_g.champion_x)
else:
print("x", end ="")
print("\nBest mass is: ", np.max(masses))
print("Worst mass is: ", np.min(masses))
best_idx = np.argmax(masses)
.......x..
Best mass is: 1344.393538122888
Worst mass is: 1300.026710328152
And we plot the trajectory found:
ax = udp_g.plot(xs[best_idx], show_gridpoints=True)
ax.view_init(90, 0)
Graident vs no gradient#
We here take a look at the difference in solving the problem using the analytical gradient and not using it. To do so we solve the same problem starting from 100 different initial guesses using both versions and compare.
from tqdm import tqdm
cpu_nog = []
cpu_g = []
fail_g = 0
fail_nog = 0
prob_g = pg.problem(udp_g)
prob_nog = pg.problem(udp_nog)
prob_g.c_tol = 1e-6
prob_nog.c_tol = 1e-6
for i in tqdm(range(100)):
pop_g = pg.population(prob_g, 1)
pop_nog = pg.population(prob_nog)
pop_nog.push_back(pop_g.champion_x)
start = time.time()
pop_g = algo.evolve(pop_g)
end = time.time()
if not prob_g.feasibility_f(pop_g.champion_f):
fail_g += 1
else: # We only record the time for successfull attempts
cpu_g.append(end - start)
start = time.time()
pop_nog = algo.evolve(pop_nog)
end = time.time()
if not prob_nog.feasibility_f(pop_nog.champion_f):
fail_nog += 1
else: # We only record the time for successfull attempts
cpu_nog.append(end - start)
print(f"Gradient: {np.median(cpu_g):.4e}s")
print(f"No Gradient: {np.median(cpu_nog):.4e}s")
print(f"\nGradient (n.fails): {fail_g}")
print(f"No Gradient (n.fails): {fail_nog}")
100%|██████████| 100/100 [02:09<00:00, 1.29s/it]
Gradient: 4.8809e-01s
No Gradient: 5.5510e-01s
Gradient (n.fails): 7
No Gradient (n.fails): 64
cpu_g = np.array(cpu_g)
cpu_nog = np.array(cpu_nog)
plt.hist(cpu_g[cpu_g < 2.2], bins=20, label="gradient", density=True, alpha=0.5)
plt.hist(cpu_nog[cpu_nog < 1.2], bins=20, label="no_gradient", density=True, alpha=0.5)
# plt.xlim([0,0.15])
plt.legend()
plt.title("point2point, 4 seg case")
plt.xlabel("CPU time")
plt.ylabel("N. occurrences")
Text(0, 0.5, 'N. occurrences')
