The Sims-Flanagan high-fidelity trajectory leg

The Sims-Flanagan trajectory leg [SF97] is implemented in pykep in the class pykep.leg.sims_flanagan_hf. The leg can be used to describe a low-thrust leg with low-fidelity as it assumes Keplerian dynamics and approximates the continuous thrust via a sequence of continuous, constant thrust arcs (zero-hold). The leg is defined by a starting position \(\mathbf x_s = [\mathbf r_s, \mathbf v_s, m_s]\), an arrival position \(\mathbf x_f = [\mathbf r_f, \mathbf v_f, m_f]\) and a time of flight \(T\).

A sequence of throttles \(\mathbf u = [u_{x0}, u_{y0}, u_{z0}, u_{x1}, u_{y1}, u_{z1}, u_{x2}, u_{y2}, u_{z2}, ... ]\) define the direction and magnitude of the continuous throttle vector along each segment (i.e. trajectory parts of equal temporal length \(\frac Tn\)).

In this tutorial we show the basic API to interface with the class pykep.leg.sims_flanagan_hf efficiently.

We start with some imports:

import pykep as pk
import numpy as np
import time
import pygmo as pg
import pygmo_plugins_nonfree as ppnf


from matplotlib import pyplot as plt
from mpl_toolkits import mplot3d

%matplotlib inline

pk.ta.get_pc(1e-16, pk.optimality_type.MASS)

C++ datatype            : double
Tolerance               : 1e-16
High accuracy           : false
Compact mode            : false
Taylor order            : 20
Dimension               : 14
Time                    : 0
State                   : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Parameters              : [0, 0, 0, 0, 0]

We then define the spacecraft propulsion system and the initial and final state. In this case they are not related to any orbital mechanics and are chosen arbitrarily for the purpose of demostrating the API.

# Problem data
mu = pk.MU_SUN
max_thrust = 0.12
isp = 3000

# Initial state
ms = 1500.0
rs = np.array([1, 0.1, -0.1]) * pk.AU
vs = np.array([0.2, 1, -0.2]) * pk.EARTH_VELOCITY

# Final state
mf = 1300.0
rf = np.array([-1.2, -0.1, 0.1]) * pk.AU
vf = np.array([0.2, -1.023, 0.44]) * pk.EARTH_VELOCITY

# Throttles and tof
nseg = 10
cut = 0.6
throttles = np.random.uniform(-1, 1, size=(nseg * 3))
tof = 324.0 * pk.DAY2SEC

Now we instantiate the leg:

# We are now ready to instantiate a leg
sf = pk.leg.sims_flanagan_hf(
    rvs=[rs, vs],
    ms=ms,
    throttles=throttles,
    rvf=[rf, vf],
    mf=mf,
    tof=tof,
    max_thrust=max_thrust,
    veff=isp * pk.G0,
    mu=mu,
    cut=cut,
)

And plot the trajectory represented by the random sequence of throttles.

# Making the axis
ax = pk.plot.make_3Daxis(figsize=(7, 7))

# Adding the Sun Earth and the boundary states
udpla = pk.udpla.jpl_lp(body="EARTH")
earth = pk.planet(udpla)
pk.plot.add_sun(ax, s=40)
pk.plot.add_planet_orbit(ax, earth)
ax.scatter(rs[0] / pk.AU, rs[1] / pk.AU, rs[2] / pk.AU, c="k", s=20)
ax.scatter(rf[0] / pk.AU, rf[1] / pk.AU, rf[2] / pk.AU, c="k", s=20)


# Plotting the trajctory leg
ax = pk.plot.add_sf_hf_leg(
    ax, sf, units=pk.AU, N=10, show_gridpoints=True, show_throttles=True, length=0.1, arrow_length_ratio=0.5
)


# Making the axis nicer
D=1
ax.set_xlim(-D,D)
ax.set_ylim(-D/2,D*3/2)
ax.view_init(90,0)
ax.axis('off');

The gradients

The main function of this leg is to serve as a fwd-bck shooting blosk in larger optimization problems. The efficent autodiff computation of the gradients of the mismatch and throttle constraints are thus provided for convenience. Indicating the mismatch constraints with \(\mathbf{mc}\) and the throttle constraints with \(\mathbf{tc}\) we seek the following quntities:

\[ \frac{\partial \mathbf {mc}}{\partial \mathbf x_s} \rightarrow (7\times7) \]

\[ \frac{\partial \mathbf {mc}}{\partial \mathbf x_f} \rightarrow (7\times7) \]

\[ \frac{\partial \mathbf {mc}}{\partial \mathbf {\tilde u}} \rightarrow (7\times(3\mathbf{nseg} + 1)) \]

\[ \frac{\partial \mathbf {tc}}{\partial \mathbf u} \rightarrow (\mathbf{nseg} \times3\mathbf{nseg}) \]

where \(\mathbf x_s=[\mathbf r, \mathbf v, m]\) is the starting state, \(\mathbf x_f\) is the final state \(\mathbf u=[u_0, u_1 ... u_{nseg}]\) are the throttles and \(\mathbf {\tilde u}=[\mathbf u, tof]\) the trottles augmented with \(tof\) .

dmcdxs, dmcdxf, dmcdu = sf.compute_mc_grad()
print("dmcdxs shape: ", dmcdu.shape)
print("dmcdxf shape: ", dmcdu.shape)
print("dmcdu shape: ", dmcdu.shape)
dtcdu = sf.compute_tc_grad()
print("dtcdu shape: ", dtcdu.shape)

dmcdxs shape:  (7, 31)
dmcdxf shape:  (7, 31)
dmcdu shape:  (7, 31)
dtcdu shape:  (10, 30)

Gradient correctness

We check the gradient correctness against numerical differentiation.

First the mismatch constraints

# These functions are to be used in an estimate_gradient_h call
def mc_from_rvms(sf, rvms):
    sf.rvms = rvms
    return sf.compute_mismatch_constraints()

def mc_from_rvmf(sf, rvmf):
    sf.rvmf = rvmf
    return sf.compute_mismatch_constraints()

def mc_from_utilda(sf, utilda):
    sf.throttles = utilda[:-1]
    sf.tof = utilda[-1]
    return sf.compute_mismatch_constraints()

# Reference values where to compute the numerical gradients
rvms = sf.rvms.copy()
rvmf = sf.rvmf.copy()
utilda = sf.throttles.copy() + [sf.tof]
# Numerical gradients
dmcdxs_numerical = pg.estimate_gradient_h(lambda x: mc_from_rvms(sf, x), rvms, 1e-6).reshape(7,7)
dmcdxf_numerical = pg.estimate_gradient_h(lambda x: mc_from_rvmf(sf, x), rvmf, 1e-6).reshape(7,7)
dmcdu_numerical = pg.estimate_gradient_h(lambda x: mc_from_utilda(sf, x), utilda, 1e-6).reshape(7, nseg*3 + 1)

We check against the analytical gradients that all went well.

print("Check gradient of mc wrt xs: ", np.allclose(dmcdxs_numerical, dmcdxs, rtol=1e-5, atol=1e-5))  # True
print("Check gradient of mc wrt xf: ", np.allclose(dmcdxf_numerical, dmcdxf, rtol=1e-5, atol=1e-5))  # True
print("Check gradient of mc wrt u:  ", np.allclose(dmcdu_numerical, dmcdu, rtol=1e-3, atol=1e-3))  # True

Check gradient of mc wrt xs:  True
Check gradient of mc wrt xf:  True
Check gradient of mc wrt u:   True

Then the throttle constraints

def tc_from_u(sf, u):
    sf.throttles = u
    return sf.compute_throttle_constraints()

# Reference values where to compute the numerical gradients
u = sf.throttles.copy()
# Numerical gradient
dtcdu_numerical = pg.estimate_gradient_h(lambda x: tc_from_u(sf, x), u, 1e-6).reshape(nseg, nseg*3)

print("Check gradient of tc wrt u:  ", np.allclose(dtcdu_numerical, dtcdu, rtol=1e-5, atol=1e-5))  # True

Check gradient of tc wrt u:   True

Custom dynamics

We can also build a pykep.leg.sims_flanagan_hf leg using a custom (zero-hold) dynamics to propagate the various segments. The dynamics needs to be written taking care of a number of things:

• It must have, as parameters, \(\mu, v_{eff}, T_1, T_2, T_3\), where these last three values represent the Cartesian component of a zero-hold thrust in some reference frame.

• It needs to have 7 dimensions with the last variable being the mass.

• Its variational version must have the full state and \(T_1, T_2, T_3\) as variational paramters.

In pykep all Taylor adaptive integrators that are called “zero_hold” respect such an interface and can thus be used.

Lets see how to use the CR3BP zero hold dynamics for this purpose.

# We are going to use a json file which contains data of problems connecting 
# periodic orbits in the cr3bp
import json
with open('sims_flanagan_hf_leg_problems.json', 'r') as f:
    problems = json.load(f)

# We create the zero-hold integrators
ta = pk.ta.get_zero_hold_cr3bp(1e-12)
ta_var = pk.ta.get_zero_hold_cr3bp_var(1e-6)
tas = [ta, ta_var]

# Change the following line to use a different problem from the json file.
problem_id = "instances_ys_mss ID2"
# We load the values in the json dictionary
state_s = problems[problem_id]["state_s"]
period_s = problems[problem_id]["period_s"]
state_f = problems[problem_id]["state_f"]
period_f = problems[problem_id]["period_f"]
max_thrust = problems[problem_id]["max_thrust"]
veff = problems[problem_id]["veff"]
mu_cr3bp = problems[problem_id]["mu_cr3bp"]

# We create new random throttles
nseg=20
throttles = np.random.uniform(-1, 1, size=(nseg * 3))

max_thrust = 4e-2

# We are now ready to instantiate a generalized leg
sf = pk.leg.sims_flanagan_hf(
    state_s=state_s,
    throttles=throttles,
    state_f=state_f,
    tof=3.141592653589793,
    max_thrust=max_thrust,
    veff=veff,
    mu=mu_cr3bp,
    cut=cut,
    tas=tas,
)

plotting is different as we are in the CR3BP!

def plot_leg_cr3bp(sf, mu_cr3bp):
    # Making the axis
    ax = pk.plot.make_3Daxis(figsize=(7, 7))

    # Adding the Sun Earth and the boundary states
    #ax.scatter(-mu_cr3bp, 0, 0)
    ax.scatter(1-mu_cr3bp, 0, 0)

    # Adding an arrow to point towards the secondary
    # Coordinates of the arrow start point
    start = [1-mu_cr3bp, 0, 0]
    # Direction vector pointing towards the secondary
    direction = [-1, 0, 0]
    # Add the arrow to the 3D plot
    ax.quiver(start[0], start[1], start[2],
            direction[0], direction[1], direction[2],
            length=0.05, color='r', normalize=True)

    # Adding the starting and end orbits
    ta.pars[:] = [mu_cr3bp, 1., 0., 0., 0.]
    # START
    ta.time = 0
    ta.state[:] = state_s
    res = ta.propagate_grid(np.linspace(0, period_s, 2000))[-1]
    ax.plot(res[:,0], res[:,1],res[:,2])
    #END
    ta.time = 0
    ta.state[:] = state_f
    res = ta.propagate_grid(np.linspace(0, period_f, 2000))[-1]
    ax.plot(res[:,0], res[:,1],res[:,2])

    # Plotting the trajctory leg
    ax = pk.plot.add_sf_hf_leg(
        ax, sf, units=1., N=10, show_gridpoints=False, show_throttles=True, length=0.01, arrow_length_ratio=0.1,
    )

    # Making the axis nicer
    ax.view_init(-30,-90)

plot_leg_cr3bp(sf, mu_cr3bp)

Embedding the leg into an optimization problem

We are now going to use the leg we instantiated above in an optimization problem (a pagmo UDP) looking for a sequence of throttles and a time of flight able to “close” the trajectory. Our decision vector will thus be:

\[ \mathbf x = [m_f, [u_x, u_y, u_z]\times n_{seg}, tof] \]

We start by writing a minimal UPD (constructor, bounds fitness). Unfortunately we need quite some lines for this as to comply with pagmo UDP requirements, but the heart of it is rather simple. Its a leg and we optimize for tof and throttles as to close the mismatches.

class p2p_in_cr3bp:
    # The Constructor -----------------------------------------------------------------------------------
    def __init__(
        self,
        state_s,
        state_f,
        nseg,
        max_thrust,
        veff,
        mu,
        cut,
        mf_bounds=[0.5, 1.0],
        tof_bounds=[1.0, 5.0],
    ):
        # We create the zero-hold integrators
        ta = pk.ta.get_zero_hold_cr3bp(1e-12)
        ta_var = pk.ta.get_zero_hold_cr3bp_var(1e-6)
        tas = [ta, ta_var]
        # We store necessary data
        self.nseg = nseg
        self.mf_bounds = mf_bounds
        self.tof_bounds = tof_bounds
        # We instantiate the leg (tof, final mass and throttles will be changed when the leg is set from the decision vector,
        # hence they are just dummy values to construct the object)
        self.leg = pk.leg.sims_flanagan_hf(
            state_s=state_s,
            throttles=[0.0] * nseg * 3,
            state_f=state_f,
            tof=2.0,
            max_thrust=max_thrust,
            veff=veff,
            mu=mu,
            cut=cut,
            tas=tas,
        )

    # The mandatory get_bounds UDP method ----------------------------------------------------------------
    def get_bounds(self):
        lb = [self.mf_bounds[0]] + [-1, -1, -1] * self.nseg + [self.tof_bounds[0]]
        ub = [self.mf_bounds[1]] + [1, 1, 1] * self.nseg + [self.tof_bounds[1]]
        return (lb, ub)

    def get_nec(self):
        return 7

    def get_nic(self):
        return self.nseg

    # Convenience method to set the leg from the decision vector
    def _set_leg_from_x(self, x):
        # x = [mf, throttles, tof]
        # We set the data in the leg using the decision vector
        self.leg.tof = x[-1]
        self.leg.mf = x[0]
        self.leg.throttles = x[1:-1]  # tof excluded

    # Fitness function -----------------------------------------------------------------------------------
    def fitness(self, x):
        # 1 - We set the leg using data in the decision vector
        self._set_leg_from_x(x)

        obj = -x[0]
        # 2 - We compute the constraints violations (mismatch+throttle)
        ceq = self.leg.compute_mismatch_constraints()
        cineq = self.leg.compute_throttle_constraints()
        retval = np.array([obj] + ceq + cineq)  # here we can sum lists
        return retval

    def plot(self, x):
        # Setting the leg
        self._set_leg_from_x(x)
        plot_leg_cr3bp(self.leg, mu_cr3bp)

We can now instatiate the UDP and hence the pagmo problem:

# Problem
# IMPORTANT: according to the problem selected the tof_bounds should be adjusted (e.g.[ 1,5] or [5,10])
tof_bounds=[15.0, 15.0]
udp = p2p_in_cr3bp(state_s, state_f, 20, max_thrust, veff, mu_cr3bp, cut = 0.6,  tof_bounds=tof_bounds)
prob = pg.problem(udp)
prob.c_tol = 1e-6 # ensures pagmo consider feasible constraint at this level (this is NOT changing SNOPT behaviour)

# Algorithm
#snopt72 = "/usr/local/lib/libsnopt7_c.so"
snopt72 = "/Users/dario.izzo/opt/libsnopt7_c.dylib"
# Activate screen output here to see what the SQP solver is doing
uda = ppnf.snopt7(library=snopt72, minor_version=2, screen_output=False)
uda.set_integer_option("Major iterations limit", 3500)
uda.set_integer_option("Minor iterations limit", 20000)
uda.set_numeric_option("Major optimality tolerance", 1e-3)
uda.set_numeric_option("Major feasibility tolerance", 1e-7) # must be low else will not be met
algo = pg.algorithm(uda)

… and call the solver (call a few times to get multiple local minima)

pop = pg.population(prob, 1)
pop = algo.evolve(pop)

udp.plot(pop.champion_x)
print("Final mass: ", -pop.champion_f[0])
print("Transfer time: ", pop.champion_x[-1])
print("Feasibility: ", np.linalg.norm(pop.champion_f[1:8]))

Final mass:  0.702380219162
Transfer time:  15.0
Feasibility:  0.574999456059

Adding the gradients

The UDP we wrote above does not have any gradient information, as one can verify printing on screen the pagmo problem object. This means that solvers, when called to solve that problem, will not have access to the gradient and must compute it numerically if they need it (SQP methods such as SNOPT do ned it).

The {class}`pykep.leg.sims_flanagan_hf, as explained above, has the possibility to return various gradients which can be used to assemble the final gradients of the problem fitness. Lets see how

def _gradient(self, x):
    self._set_leg_from_x(x)
    _, mcg_xf, mcg_th_tof = self.leg.compute_mc_grad()
    tcg_th = self.leg.compute_tc_grad()

    # 1 - The gradient of the objective function (obj = -mf)
    retval = [-1.0]
    # 2 - The gradient of the mismatch contraints (mcg). We divide them in pos, vel mass as they have different scaling units
    # pos
    for i in range(7):
        # First w.r.t. mf
        retval.append(mcg_xf[i, -1])
        # Then the [throttles, tof]
        retval.extend(mcg_th_tof[i, :])
    # 3 -  The gradient of the throttle constraints
    for i in range(self.nseg):
        retval.extend(tcg_th[i, 3 * i : 3 * i + 3])
    return retval

def _gradient_sparsity(self):
    dim = 2 + 3 * self.nseg
    # The objective function only depends on the final mass, which is the first variable in the decision vector.
    retval = [[0, 0]]
    # The mismatch constraints depend on all variables.
    for i in range(1, 8):
        for j in range(dim):
            retval.append([i, j])
    # The throttle constraints only depend on the specific throttles (3).
    for i in range(self.nseg):
        retval.append([8 + i, 3 * i + 1])
        retval.append([8 + i, 3 * i + 2])
        retval.append([8 + i, 3 * i + 3])
    # We return the sparsity pattern
    return retval

#Monkey patch the class
p2p_in_cr3bp.gradient = _gradient
p2p_in_cr3bp.gradient_sparsity = _gradient_sparsity

# Problem
udp = p2p_in_cr3bp(state_s, state_f, 21, max_thrust, veff, mu_cr3bp, cut = 0.6,  tof_bounds=tof_bounds)
prob = pg.problem(udp)
prob.c_tol = 1e-6

pop.champion_f[0]

-0.55276328450279788

for i in range(50):
    pop = pg.population(prob, 1)
    pop = algo.evolve(pop)
    np.linalg.norm(pop.champion_f[1:7])
    if  np.linalg.norm(pop.champion_f[1:7]) < 1e-6:
        print(f"{pop.champion_x[0]}, {np.linalg.norm(pop.champion_f[1:7])}")
        if pop.champion_f[0] < - 0.97:
            break

0.9577693376520025, 2.7591078790100513e-07
0.9522004564482348, 4.3857911616309804e-08
0.9523170432472177, 1.2203780182075033e-08
0.9669911466576542, 1.3130708277853673e-08
0.9541984845512791, 4.949575519596202e-07
0.9522698624434427, 6.180347676184902e-08
0.9523341681149741, 2.8965476229862366e-08
0.9523380938970615, 5.4430773066174426e-08
0.952304012313531, 3.5039231854616343e-09
0.6463524868386623, 6.0732262104787515e-09
0.9522194438942301, 1.0034133455068918e-08
0.95041198129228, 2.6799413340676307e-10
0.9570591541651174, 1.6048735189261368e-09
0.9523379340077337, 5.344843320803091e-08
0.9672314455088733, 5.508248571390954e-07
0.952337860926775, 2.7567865283998857e-07
0.6463525119756856, 1.592713967428314e-07
0.9515285857622272, 3.914017358962119e-08
0.9523380446079853, 6.078633909452506e-09
0.9522849129952831, 1.6800706801766615e-07
0.948726077288913, 6.031807842010098e-12
0.952330070448648, 1.555536459649651e-07
0.9522912193286862, 1.3787207900912312e-08
0.9585120919268082, 1.0462398184310548e-07
0.9837407082269952, 2.4959347602709045e-07

udp.plot(pop.champion_x)
print("Final mass: ", -pop.champion_f[0])
print("Transfer time: ", pop.champion_x[-1])
print("Feasibility: ", np.linalg.norm(pop.champion_f[1:8]))

Final mass:  0.983740708227
Transfer time:  15.0
Feasibility:  2.51763032299e-07

max_thrust

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 max_thrust

NameError: name 'max_thrust' is not defined