Trajectory Optimization

In pykep both direct and indirect optimization transcriptions as well as evolutionary encodings are provided to perform spacecraft trajectory optimization. Most direct techniques provided are variations to the “Sims” transcription [SF97], while indirect techniques are based on some version of Pontryagin Maximum Principle. The evolutionary encoding are mostly based on the work performed at ESA’ Advanced Concepts Team ([Izz10], [ISM+13]).

Most of the classes in this Trajectory Optimization module are provided as optimization problems (NLPs) in the for of User Defined Problems (UDP) compatible with the pygmo [BI20] python package.

Direct

When solving optimal control problems (OCPs), the unknown is a function (i.e. the control), hence one is dealing with an infinite dimensional search space. A common approach to overcome this difficulty is to use direct methods. These are based on the idea of discretizing the problem introducing some time grid and thus transforming the OCP into a NLP (Non Linear Programming) problem to then solve it using numerical solvers available for the task.

Sims-Flanagan

One of pykep’s historical favourite transcription methods is the Sims [SF97] transcription which is called Sims-Flanagan as per the authorship of the first brief conferecne paper introducing it (anecdotically the work was, though, entirely done by Sims).

---

class pykep.trajopt.sf_point2point(rvs=[_np.array([1, 0.1, -0.1]) * pk.AU, _np.array([0.2, 1, -0.2]) * pk.EARTH_VELOCITY], rvf=[_np.array([-1.2, -0.1, 0.1]) * pk.AU, _np.array([0.2, -1.023, 0.44]) * pk.EARTH_VELOCITY], ms=1000., mu=pk.MU_SUN, max_thrust=0.12, veff=3000 * pk.G0, tof_bounds=[80., 400.], mf_bounds=[200.0, 1000.0], nseg=10, cut=0.6, mass_scaling=1000, r_scaling=pk.AU, v_scaling=pk.EARTH_VELOCITY, with_gradient=True)

Represents the optimal low-thrust transfer between two fixed points using the Sims-Flanagan (direct) method.

This problem works internally using the sims_flanagan and manipulates its transfer time T, final mass mf and the controls as to link two fixed points in space with a low-thrust trajectory.

It can be used to better profile and understand performances of optimizers on this type of direct approach, but has a limited use in the design of interplanetary trajectories as per the fixed point limitation.

The decision vector is:

z = [mf, throttles, tof]

where throttles is a vector of throttles structures as [u0x, u0y,u0z, …]. By throttles we intend non dimensiona thrust levels in [0,1].

Initializes the sf_point2point instance with given parameters.

Args:

rvs (list): Initial position and velocity vectors. Defaults to two vectors scaled by AU and Earth’s velocity.

rvf (list): Final position and velocity vectors. Defaults to two vectors scaled by AU and Earth’s velocity.

ms (float): Initial spacecraft mass in kg. Defaults to 1000 kg.

mu (float): Gravitational parameter, default is for the Sun (MU_SUN).

max_thrust (float): Maximum thrust in Newtons. Defaults to 0.12 N.

isp (float): Specific impulse in seconds. Defaults to 3000 s.

tof_bounds (list): Bounds for time of flight in days. Defaults to [0, 400] days.

mf_bounds (list): Bounds for final mass in kg. Defaults to [200.0, 1000.0] kg.

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

mass_scaling (float): Scaling factor for mass (used to scale constraints). Defaults to 1000.

r_scaling (float): Scaling factor for distance, (used to scale constraints). Defaults AU (AU).

v_scaling (float): Scaling factor for velocity (used to scale constraints). Defaults the Earth’s velocity (EARTH_VELOCITY).

with_gradient (bool): Indicates if gradient information should be used. Defaults True.

plot(x, ax=None, units=pk.AU, show_midpoints=False, show_gridpoints=False, show_throttles=False, length=0.1, arrow_length_ratio=0.05, **kwargs)

Plots the trajectory leg 3D axes.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to _pk.AU.

show_midpoints (bool, optional): Whether to show midpoints on the trajectory. Defaults to False.

show_gridpoints (bool, optional): Whether to show grid points on the trajectory. Defaults to False.

show_throttles (bool, optional): Whether to show throttle vectors. Defaults to False.

length (float, optional): Length of the throttle vectors. Defaults to 0.1.

arrow_length_ratio (float, optional): Arrow length ratio for the throttle vectors. Defaults to 0.05.

**kwargs: Additional keyword arguments for the plot.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis with the plotted trajectory.

pretty(x)

Prints a detailed representation of the Point to point problem.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

---

class pykep.trajopt.sf_pl2pl(pls=pk.planet(pk.udpla.jpl_lp(body='EARTH')), plf=pk.planet(pk.udpla.jpl_lp(body='MARS')), ms=1500, mu=pk.MU_SUN, max_thrust=0.12, veff=3000 * pk.G0, t0_bounds=[6700.0, 6800.0], tof_bounds=[200.0, 300.0], mf_bounds=[1300.0, 1500.0], vinfs=3.0, vinff=0.0, nseg=10, cut=0.6, mass_scaling=1500, r_scaling=pk.AU, v_scaling=pk.EARTH_VELOCITY, with_gradient=True)

Represents the optimal low-thrust transfer between two planet using the Sims-Flanagan (direct) method.

This problem works internally using the sims_flanagan and manipulates its initial and final states, as well as its transfer time T, final mass mf and the controls as to link the two planets with a low-thrust trajectory.

The particular transcription used is suitable only for few revolutions, after which convergence will start to be problematic.

The decision vector is:

z = [t0, mf, Vsx, Vsy, Vsz, Vfx, Vfy, Vfz, throttles, tof] - all in S.I. units except t0 and tof in days

where throttles is a vector of throttles structured as [u0x, u0y,u0z, …]. By throttles we intend non dimensiona thrust levels in [0,1].

Args:

pls (planet): Initial planet. Defaults to jpl_lp Earth.

plf (planet): Final planet. Defaults to jpl_lp Mars.

ms (float): Initial spacecraft mass in kg. Defaults to 1000 kg.

mu (float): Gravitational parameter, default is for the Sun (MU_SUN).

max_thrust (float): Maximum thrust in Newtons. Defaults to 0.12 N.

veff (float): Exhaust velocity (Specific impulse * _pk.G0) in m/s. Defaults to 3000 s * _pk.G0.

t0_bounds (list): Bounds for departure epoch in MJD2000. Defaults to [6700.0, 6800.0].

tof_bounds (list): Bounds for time of flight in days. Defaults to [200, 300] days.

mf_bounds (list): Bounds for final mass in kg. Defaults to [1300.0, 1500.0] kg.

vinfs (float): Allowed magnitude for the departure’s relative velocity in km/s. Defaults to 3.

vinff (float): Allowed magnitude for the arrival’s relative velocity in km/s. Defaults to 0.

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

mass_scaling (float): Scaling factor for mass (used to scale constraints). Defaults to 1500.

r_scaling (float): Scaling factor for distance, (used to scale constraints). Defaults AU (AU).

v_scaling (float): Scaling factor for velocity (used to scale constraints). Defaults the Earth’s velocity (EARTH_VELOCITY).

with_gradient (bool): Indicates if gradient information should be used. Defaults True.

plot(x, ax=None, units=pk.AU, show_midpoints=False, show_gridpoints=False, show_throttles=False, length=0.1, arrow_length_ratio=0.05, **kwargs)

Plots the trajectory leg 3D axes.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to _pk.AU.

show_midpoints (bool, optional): Whether to show midpoints on the trajectory. Defaults to False.

show_gridpoints (bool, optional): Whether to show grid points on the trajectory. Defaults to False.

show_throttles (bool, optional): Whether to show throttle vectors. Defaults to False.

length (float, optional): Length of the throttle vectors. Defaults to 0.1.

arrow_length_ratio (float, optional): Arrow length ratio for the throttle vectors. Defaults to 0.05.

**kwargs: Additional keyword arguments for the plot.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis with the plotted trajectory.

pretty(x)

Prints a human readable representation of the transfer.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

---

class pykep.trajopt.sf_pl2pl_alpha(pls=pk.planet(pk.udpla.jpl_lp(body='EARTH')), plf=pk.planet(pk.udpla.jpl_lp(body='MARS')), ms=1500, mu=pk.MU_SUN, max_thrust=0.12, veff=3000 * pk.G0, t0_bounds=[6700.0, 6800.0], tof_bounds=[200.0, 300.0], mf_bounds=[1300.0, 1500.0], vinfs=3.0, vinff=0.0, nseg=10, cut=0.6, mass_scaling=1500, r_scaling=pk.AU, v_scaling=pk.EARTH_VELOCITY)

Represents the optimal low-thrust transfer between two planet using the Sims-Flanagan (direct) method.

This problem works internally using the sims_flanagan and manipulates its initial and final states, as well as its transfer time T, final mass mf and the controls as to link the two planets with a low-thrust trajectory.

The particular transcription used is suitable only for few revolutions, after which convergence will start to be problematic.

The decision vector is:

z = [t0, mf, Vsx, Vsy, Vsz, Vfx, Vfy, Vfz, alphas, throttles, tof] - all in S.I. units except t0 and tof in days

where throttles is a vector of throttles structured as [u0x, u0y,u0z, …]. By throttles we intend non dimensiona thrust levels in [0,1].

Args:

pls (planet): Initial planet. Defaults to jpl_lp Earth.

plf (planet): Final planet. Defaults to jpl_lp Mars.

ms (float): Initial spacecraft mass in kg. Defaults to 1000 kg.

mu (float): Gravitational parameter, default is for the Sun (MU_SUN).

max_thrust (float): Maximum thrust in Newtons. Defaults to 0.12 N.

veff (float): Exhaust velocity (Specific impulse * _pk.G0) in m/s. Defaults to 3000 s * _pk.G0.

t0_bounds (list): Bounds for departure epoch in MJD2000. Defaults to [6700.0, 6800.0].

tof_bounds (list): Bounds for time of flight in days. Defaults to [200, 300] days.

mf_bounds (list): Bounds for final mass in kg. Defaults to [1300.0, 1500.0] kg.

vinfs (float): Allowed magnitude for the departure’s relative velocity in km/s. Defaults to 3.

vinff (float): Allowed magnitude for the arrival’s relative velocity in km/s. Defaults to 0.

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

mass_scaling (float): Scaling factor for mass (used to scale constraints). Defaults to 1500.

r_scaling (float): Scaling factor for distance, (used to scale constraints). Defaults AU (AU).

v_scaling (float): Scaling factor for velocity (used to scale constraints). Defaults the Earth’s velocity (EARTH_VELOCITY).

plot(x, ax=None, units=pk.AU, show_midpoints=False, show_gridpoints=False, show_throttles=False, length=0.1, arrow_length_ratio=0.05, **kwargs)

Plots the trajectory leg 3D axes.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to _pk.AU.

show_midpoints (bool, optional): Whether to show midpoints on the trajectory. Defaults to False.

show_gridpoints (bool, optional): Whether to show grid points on the trajectory. Defaults to False.

show_throttles (bool, optional): Whether to show throttle vectors. Defaults to False.

length (float, optional): Length of the throttle vectors. Defaults to 0.1.

arrow_length_ratio (float, optional): Arrow length ratio for the throttle vectors. Defaults to 0.05.

**kwargs: Additional keyword arguments for the plot.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis with the plotted trajectory.

pretty(x)

Prints a human readable representation of the transfer.

Args:

x (list): The decision vector containing final mass, thrust direction, and time of flight.

Zero-Order-Hold

The pykep original extension to the original Sims-Flanagan model is called the Zero-Order_hold transcription mehtod as it substitutes impulses with constant control segments. It has been developed and constantly refined throughout years starting from the early appearence in 2010 in the work from Yam et al. [YIB10] who reported it as a high-fidelity version of the original transcription.

class pykep.trajopt.zoh_point2point(states=[1.2, 0.0, -0.01, 0.01, 1.0, -0.01], statef=[1.0, 0.0, -0.0, 0.01, 1.1, -0.0], ms=1.0, max_thrust=0.22, tof_bounds=[3.4, 8.6], mf_bounds=[0.2, 1], nseg=10, cut=0.6, tas=(_pk.ta.get_zoh_kep(1e-10), None), time_encoding='uniform', w_bounds_softmax=[-1.0, 1.0])

Represents the optimal low-thrust transfer between two fixed points using pykep’s Zero Order Hold (direct) trajectory legs.

This problem works internally using the sims_flanagan and manipulates its transfer time T, final mass mf and the controls as to link two fixed points in space with a low-thrust trajectory.

It can be used to better profile and understand performances of optimizers on this type of direct approach, but has a limited use in the design of interplanetary trajectories as per the fixed point limitation.

The decision vector is:

x = [mf] + controls + tof (+ [weights for softmax])

where controls is a vector of control parameters \([T, i_x, i_y, i_z] \times n_\text{seg}\) representing magnitude and direction of the thrust applied in each segment.

Initializes the zoh_point2point instance with given parameters.

Args:

states (list): Initial state (only the first 6 states, i.e. no mass). Units as expected by the numerical integrator, defaults to [1.2, 0.0, -0.01, 0.01, 1.0, -0.01].

statef (list): Final state (only the first 6 states, i.e. no mass). Units as expected by the numerical integrator, defaults to [1.0, 0.0, -0.0, 0.01, 1.1, -0.0].

ms (float): Initial mass. Units as expected by the numerical integrator, defaults to 1.

max_thrust (float): Maximum thrust. Units as expected by the numerical integrator, defaults to 0.22.

tof_bounds (list): Bounds for time of flight. Units as expected by the numerical integrator, defaults to [3.4, 8.6].

mf_bounds (list): Bounds for final mass. Units as expected by the numerical integrator, defaults to [0.2, 1].

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

tas (tuple): (ta, ta_var) Taylor-adaptive integrators

• ta: Nominal dynamics (state dim 7, pars ≥ 4)

• ta_var: Variational dynamics (state dim 84, same pars). When None, no gradients will be used.

plot(x, ax=None, N=30, to_cartesian=lambda state: ..., mark_segments=True, mark_mismatch=True, **kwargs)

Plots the trajectory of the zero order hold point to point problem.

Args:

x (list): The decision vector containins: final mass, thrust direction, time of flight and (if time encoding is softmax) the weights for the softmax time grid.

ax (matplotlib.axes.Axes): The matplotlib axes to plot on. If None, a new figure and axes will be created.

N (int): The number of points to plot along the trajectory.

to_cartesian (Callable): A function that converts whatever state is used in the internal Taylor integrator to Cartesian (r,v).

mark_segments (bool): adds markers ath each segment edge

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The modified Axes object with the Lambert’s problem trajectory added.

pretty(x)

Prints a detailed representation of the zero order hold point to point problem.

Args:

x (list): The decision vector containins: final mass, thrust direction, time of flight and (if time encoding is softmax) the weights for the softmax time grid.

---

class pykep.trajopt.zoh_pl2pl(pls=_pk.planet(_pk.udpla.jpl_lp(body='EARTH')), plf=_pk.planet(_pk.udpla.jpl_lp(body='MARS')), ms=1.0, max_thrust=0.22, t0_bounds=[6700.0, 6800.0], tof_bounds=[3.4, 8.6], mf_bounds=[0.2, 1.0], vinf_dep_bounds=[0.0, 0.2], vinf_arr_bounds=[0.0, 0.2], nseg=10, cut=0.6, tas=(_pk.ta.get_zoh_kep(1e-10), None), time_encoding='uniform', w_bounds_softmax=[-1.0, 1.0], L=_pk.AU, V=_pk.EARTH_VELOCITY)

Represents the optimal low-thrust transfer between two planet using the Zero Order Hold (direct) method with free departure and arrival velocities.

This problem works internally using the zoh and manipulates its initial and final states, as well as its transfer time T, final mass mf and the controls as to link the two planets with a low-thrust trajectory. The spacecraft can depart and arrive with velocities different from the planets.

The problem works in non-dimensional units internally, while planets return SI units. The user must provide scaling factors (L, V, MASS) for proper conversion.

The decision vector is:

x = [t0_fractional, mf, vinf_dep, idep_x, idep_y, idep_z, vinf_arr, iarr_x, iarr_y, iarr_z] + controls + [tof] (+ [weights] for softmax)

where: - t0_fractional is in non-dimensional (fraction of the t0 bounds) - mf is non-dimensional (scaled by MASS) - vinf_dep, vinf_arr are non-dimensional (scaled by V) - idep, iarr are unit direction vectors - controls = [T, i_x, i_y, i_z] × nseg (T non-dimensional, directions unit vectors) - tof is non-dimensional (scaled by TIME = sqrt(L³/MU))

Initializes the zoh_pl2pl_free_v instance with given parameters.

Args:

pls (planet): Initial planet. Defaults to jpl_lp Earth.

plf (planet): Final planet. Defaults to jpl_lp Mars.

ms (float): Initial spacecraft mass (non-dimensional). Defaults to 1.0.

max_thrust (float): Maximum thrust (non-dimensional). Defaults to 0.22.

t0_bounds (list): Bounds for departure epoch in MJD2000. Defaults to [6700.0, 6800.0].

tof_bounds (list): Bounds for time of flight (non-dimensional). Defaults to [3.4, 8.6].

mf_bounds (list): Bounds for final mass (non-dimensional). Defaults to [0.2, 1.0].

vinf_dep_bounds (list): Bounds for departure excess velocity magnitude (non-dimensional). Defaults to [0.0, 0.2].

vinf_arr_bounds (list): Bounds for arrival excess velocity magnitude (non-dimensional). Defaults to [0.0, 0.2].

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the zoh. Defaults to 0.6.

tas (tuple): (ta, ta_var) Taylor-adaptive integrators

• ta: Nominal dynamics (state dim 7, pars ≥ 4)

• ta_var: Variational dynamics (state dim 84, same pars). When None, no gradients will be used.

time_encoding (str): Time encoding scheme. Defaults to ‘uniform’. Options: ‘uniform’, ‘softmax’.

w_bounds_softmax (list): Bounds for softmax weights. Defaults to [-1., 1.].

L (float): Length units. Defaults to AU. (user must ensure consistency with the tas assumptions)

V (float): Velocity units. Defaults to EARTH_VELOCITY. (user must ensure consistency with the tas assumptions)

plot(x, ax=None, N=30, mark_segments=True, mark_mismatch=True, **kwargs)

Plots the trajectory of the zero order hold point to point problem.

Args:

x (list): The decision vector containins: final mass, thrust direction, time of flight and (if time encoding is softmax) the weights for the softmax time grid.

ax (matplotlib.axes.Axes): The matplotlib axes to plot on. If None, a new figure and axes will be created.

N (int): The number of points to plot along the trajectory.

mark_segments (bool): adds markers ath each segment edge

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The modified Axes object with the Lambert’s problem trajectory added.

pretty(x)

Prints a detailed representation of the zero order hold planet to planet problem with free velocities.

---

class pykep.trajopt.zoh_ss_point2point(states=[1.2, 0.0, -0.01, 0.01, 1.0, -0.01], statef=[1.0, 0.0, -0.0, 0.01, 1.1, -0.0], tof_bounds=[3.4, 8.6], nseg=10, cut=0.6, tas=(_pk.ta.get_zoh_kep(1e-10), None), time_encoding='uniform', w_bounds_softmax=[-1.0, 1.0])

Represents the time optimal Solar Sail transfer between two fixed points using pykep’s Zero Order Hold (direct) trajectory legs.

This problem works internally using the zoh_ss and manipulates its transfer time tof and the sail clock and cone angles (controls) as to link two fixed points in space with a Solar Sail trajectory.

It can be used to better profile and understand performances of optimizers on this type of direct approach, but has a limited use in the design of interplanetary trajectories as per the fixed point limitation.

The decision vector is:

x = controls + [tof] (+ [weights for softmax])

where controls is a vector of control parameters \([\alpha, \beta] \times n_\text{seg}\) representing the sail cone and clock angles, assumed fixed in the RTN frame along each segment.

Initializes the zoh_ss_point2point instance with given parameters.

Args:

states (list): Initial state. Units as expected by the numerical integrator, defaults to [1.2, 0.0, -0.01, 0.01, 1.0, -0.01].

statef (list): Final state. Units as expected by the numerical integrator, defaults to [1.0, 0.0, -0.0, 0.01, 1.1, -0.0].

tof_bounds (list): Bounds for time of flight. Units as expected by the numerical integrator, defaults to [3.4, 8.6].

nseg (int): Number of segments for the trajectory. Defaults to 10.

cut (float): Cut parameter for the sims_flanagan. Defaults to 0.6.

tas (tuple): (ta, ta_var) Taylor-adaptive integrators

• ta: Nominal dynamics (state dim 6, pars ≥ 2)

• ta_var: Variational dynamics (54, same pars). When None, no gradients will be used.

plot(x, ax=None, N=30, to_cartesian=lambda state: ..., mark_segments=True, mark_mismatch=True, plot_sail=True, sail_size=0.05, **kwargs)

Plots the trajectory of the zero order hold point to point problem.

Args:

x (list): The decision vector containing: final mass, thrust direction,

time of flight and (if time encoding is softmax) the weights for the softmax time grid.

ax (matplotlib.axes.Axes): The matplotlib axes to plot on.

If None, a new figure and axes will be created.

N (int): The number of points to plot along the trajectory.

to_cartesian (Callable): A function that converts whatever

state is used in the internal Taylor integrator to Cartesian (r,v).

mark_segments (bool): Adds markers at each segment edge.

mark_mismatch (bool): Marks the terminal mismatch point.

plot_sail (bool): Adds a visualizatio of a rectangular sail.

sail_size (float): Half-side length of the rendered sail square,

in the same units as the trajectory positions.

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The modified Axes object with the Lambert’s problem trajectory added.

pretty(x)

Prints a detailed representation of the zero order hold point to point problem.

Args:

x (list): The decision vector containins: final mass, thrust direction, time of flight and (if time encoding is softmax) the weights for the softmax time grid.

Indirect

Indirect methods are based on the Pontryagin Maximum Principle (PMP), which provides necessary conditions for optimality. These methods involve deriving the optimal control laws and the corresponding state trajectories by solving a two point boundary value problem (TPBVP) derived from the PMP.

class pykep.trajopt.pontryagin_cartesian_mass(posvel0=_posvel0, posvelf=_posvelf, tof=250, mu=pk.MU_SUN, lambda0=1.0, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=pk.AU, MU=pk.MU_SUN, MASS=1500, with_gradient=False, taylor_tolerance=1e-16, taylor_tolerance_var=1e-8)

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single-shooting indirect method for the minimum mass, fixed time, optimization of a low-thrust trajectory.

The decision vector is:

x = [lx, ly, lz, lvx, lvy, lvz, lm, l0]

Note

Gradients are also provided via the variational version of the augmented dynamics.

pykep.trajopt.pontryagin_cartesian(start=default, target=default, tof=250, mu=1.32712440018e+20, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=1.495978707e+11, TIME=31557600.0, MASS=1500, with_gradient=False)

Args:

posvel0 (list): the initial position and velocity of the spacecraft.

posvelf (list): the final position and velocity of the spacecraft.

tof (float): the time of flight of the trajectory (in days).

mu (float): the gravitational parameter of the central body.

lambda0 (float or None): multiplicative factor for the objective. If None, lambda0 is added to the decision vector and the constraint ||lambda||=1 is added so that all costates must be in [-1,1].

eps (float): the accuracy of the numerical integrator.

T_max (float): the maximum thrust of the spacecraft.

Isp (float): the specific impulse of the spacecraft.

m0 (float): the initial mass of the spacecraft.

L (float): units for length. Default is the astronomical unit (AU).

MU (float): units for the gravitational parameter of the central body. Default is the gravitational parameter of the Sun.

MASS (float): units for mass. Default is 1500 kg.

with_gradient (bool): whether to use the gradient of the constraints or not.

taylor_tolerance (float): the tolerance for the Taylor integrator.

taylor_tolerance_var (float): the tolerance for the variational Taylor integrator.

plot(x, N=100, ax3D=None)

This function plots the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot. ax3D (matplotlib.axes._axes.Axes): the axis to use for the plot. If None, a new axis is created.

Returns:

ax3D (matplotlib.axes._axes.Axes): the axis of the plot.

plot_misc(x, N=100)

This function plots the throttle, thrust direction, switching function, mass costate and Hamiltonian of the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot.

Returns:

axs (list): the list of axis of the plots.

class pykep.trajopt.pontryagin_cartesian_time(start=default, target=default, tof=250, mu=1.32712440018e+20, eps=1e-4, T_max=0.6, Isp=3000, m0=1500, L=1.495978707e+11, TIME=31557600.0, MASS=1500, with_gradient=False)

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single-shooting indirect method for the minimum time, optimization of a low-thrust trajectory.

The decision vector is:

x = [lx, ly, lz, lvx, lvy, lvz, lm, lJ, tof]

Note

Gradients are also provided via the variational version of the augmented dynamics.

Args:

source (list): the initial planet.

target (list): the final planet.

t0 (epoch): the initial epoch.

tof_guess (float): a guess for the time of flight. Bound will be defined as tof_guess/2 and 2*tof_guess.

lambda0 (float or None): multiplicative factor for the objective. If None, lambda0 is added to the decision vector and the constraint ||lambda||=1 is added so that all costates must be in [-1,1].

T_max (float): the maximum thrust of the spacecraft.

Isp (float): the specific impulse of the spacecraft.

m0 (float): the initial mass of the spacecraft.

L (float): units for length. Default is the astronomical unit (AU).

MU (float): units for the gravitational parameter of the central body. Default is the gravitational parameter of the Sun.

MASS (float): units for mass. Default is 1500 kg.

with_gradient (bool): whether to use the gradient of the constraints or not.

taylor_tolerance (float): the tolerance for the Taylor integrator.

taylor_tolerance_var (float): the tolerance for the variational Taylor integrator.

plot(x, ax3D=None, N=100, **kwargs)

This function plots the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. ax3D (matplotlib.axes._axes.Axes): the axis to use for the plot. If None, a new axis is created. N (int): the number of points to use in the plot. **kwargs: Additional keyword arguments to pass to the trajectory plot.

Returns:

ax3D (matplotlib.axes._axes.Axes): the axis of the plot.

plot_misc(x, N=100, **kwargs)

This function plots the throttle, switching function and Hamiltonian of the trajectory encoded in the decision vector x.

Args:

x (list): the decision vector. N (int): the number of points to use in the plot. **kwargs: Additional keyword arguments to pass to the plt.subplots function.

Returns:

axs (list): the list of axis of the plots.

---

Evolutionary Encodings

Some type of interplanetary trajectories can be evolved: shocking?. This AI approach to trajectory design resulted in two silver ([TAF+05], [PGJ+18]) and one gold ([ISM+13]) Humies award for human-competitive results that were produced by any form of genetic and evolutionary computation. In pykep we offer some classes that help instantiating these type of transfers amenable to evolutionary techniques.

---

class pykep.trajopt.mga(seq, t0, tof, vinf, multi_objective=False, tof_encoding="direct"", orbit_insertion=False, e_target=None, rp_target=None)

The Multiple Gravity Assist (MGA) encoding of an interplanetary trajectory.

This class may be used as a User Defined Problem (UDP) for the pygmo (http://esa.github.io/pygmo/) optimisation suite.

The decision vector (chromosome) is:

direct encoding: z = [t0, T1, T2 ... ] in [mjd2000, days, days ... ]
alpha encoding: z = [t0, T, a1, a2 ...] in [mjd2000, days, nd, nd ... ]
eta encoding: z = [t0, n1, n2, n3 ...] in [mjd2000, nd, nd ...]

Note

The time of flights of a MGA trajectory (and in general) can be encoded in different ways. When they are directly present in the decision vector, we talk about a direct encoding. This is the most ‘evolvable’ encoding but also the one that requires the most problem knowledge (e.g. to define the bounds on each leg) and is not very flexible in dealing with constraints on the total time of flight. The alpha and eta encodings, instead, allow to only specify bounds on the time of flight of the entire trajectory, and not on the single legs: a property that is attractive for multi-objective optimization, for example.

In the alpha encoding each leg time-of-flight is decoded as follows,

\(T_i = T \log\alpha_i / \sum_n\log\alpha_n\).

In the eta encoding each leg time-of-flight is decoded as follows,

\(T_i = (T_{max} - \sum_{j=0}^{(i-1)}T_j) \eta_i\)

The chromosome dimension for the direct and eta encoding is the same, while the alpha encoding requires one more gene.

Note

The resulting problem is box-bounded (unconstrained).

Args:

seq (list [planet]): sequence of planetary encounters including the departure body.

t0 (list [float or epoch]): lower and upper bounds for the launch epoch. When floats are used MJD2000 is assumed.

tof (list or float): defines the bounds on the time of flight. If tof_encoding is ‘direct’, this contains a list of 2D lists defining the upper and lower bounds on each leg. If tof_encoding is ‘alpha’, this contains a 2D list with the lower and upper bounds on the total time-of-flight. If tof_encoding is ‘eta’ tof is a float defining an upper bound for the time-of-flight.

vinf (float): the vinf provided at launch for free

multi_objective (bool): when True constructs a multiobjective problem (dv, T). In this case, ‘alpha’ or eta encodings are recommended

tof_encoding (str): one of ‘direct’, ‘alpha’ or ‘eta’. Selects the encoding for the time of flights

orbit_insertion (bool): when True the arrival dv is computed as that required to acquire a target orbit defined by e_target and rp_target

e_target (float): if orbit_insertion is True this defines the target orbit eccentricity around the final planet

rp_target (float): if orbit_insertion is True this defines the target orbit pericenter around the final planet (in m)

max_revs (int): maximal number of revolutions for lambert transfer

plot(x, ax=None, units=pk.AU, N=60, c_orbit='dimgray', c_lambert='indianred', leg_ids=[], figsize=(5, 5), **kwargs)

Plots the trajectory encoded into x in 3D axes.

Args:

x (list): The decision vector in the correct tof encoding.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to pk.AU.

N (int, optional): The number of points to use when plotting the trajectory. Defaults to 60.

c_orbit (str, optional): The color of the planet orbits. Defaults to ‘dimgray’.

c (str, optional): The color of the trajectory. Defaults to ‘indianred’.

figsize (tuple): The figure size (only used if a*ax* is None and axis have to be created.), Defaults to (5, 5).

leg_ids (list): selects the legs to plot. Optional, defaults to all legs.

**kwargs: Additional keyword arguments to pass to the trajectory plot (all Lambert arcs)

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis where the trajectory was plotted.

pretty(x)

Prints a human readable representation of the transfer.

Args:

x (list): The decision vector in the correct tof encoding.

class pykep.trajopt.mga_1dsm(seq, t0, tof, vinf_bounds=[0.5, 2.5], add_vinf_dep=False, add_vinf_arr=True, tof_encoding='direct', multi_objective=False, orbit_insertion=False, e_target=None, rp_target=None, eta_bounds=[0.01, 0.99], rp_ub=30)

This class transcribes a Multiple Gravity Assist trajectory with one deep space maneuver per leg (MGA-1DSM) into an optimisation problem. It may be used as a User Defined Problem (UDP) for the pygmo (http://esa.github.io/pygmo/) optimisation suite.

• Izzo, Dario. “Global optimization and space pruning for spacecraft trajectory design.” Spacecraft Trajectory Optimization 1 (2010): 178-200.

The decision vector (chromosome) is:

direct encoding: [t0] + [u, v, Vinf, eta1, T1] + [beta, rp/rV, eta2, T2] + ...
alpha encoding:  [t0] + [u, v, Vinf, eta1, a1] + [beta, rp/rV, eta2, a2] + ... + [T]
eta encoding:    [t0] + [u, v, Vinf, eta1, n1] + [beta, rp/rV, eta2, n2] + ...

where t0 is a mjd2000, Vinf is in km/s, T in days, beta in radians and the rest non dimensional.

Note

The time of flights of a MGA-1DSM trajectory (and in general) can be encoded in different ways. When they are directly present in the decision vector, we have a direct encoding. This is the most ‘evolvable’ encoding but also the one that requires the most problem knowledge (e.g. to define the bounds on each leg) and is not very flexible in dealing with constraints on the total time of flight. The alpha and eta encodings, instead, allow to only specify bounds on the time of flight of the entire trajectory, and not on the single legs: a property that is attractive for multi-objective optimization, for example.

In the alpha encoding each leg time-of-flight is decoded as follows, T_i = T log(alpha_i) / sum_n(log(alpha_n)). In the eta encoding each leg time-of-flight is decoded as follows, T_i = (tof_max - sum_0^(i-1)(T_j)) * eta_i

The chromosome dimension for the direct and eta encoding is the same, while the alpha encoding requires one more gene.

Note

The resulting problem is box-bounded (unconstrained).

seq (list [planet]): sequence of planetary encounters including the departure body.

t0 (list [float or epoch]): lower and upper bounds for the launch epoch. When floats are used MJD2000 is assumed.

tof (list or float): defines the bounds on the time of flight. If tof_encoding is ‘direct’, this contains a list of 2D lists defining the upper and lower bounds on each leg. If tof_encoding is ‘alpha’, this contains a 2D list with the lower and upper bounds on the total time-of-flight. If tof_encoding is ‘eta’, tof is a float defining an upper bound for the time-of-flight.

vinf_bounds (list [float]): the minimum and maximum allowed initial hyperbolic velocity (at launch), in km/sec.

add_vinf_dep (bool): when True the computed Dv includes the initial hyperbolic velocity (at launch).

add_vinf_arr (bool): when True the computed Dv includes the final hyperbolic velocity (at the last planet).

tof_encoding (str): one of ‘direct’, ‘alpha’ or ‘eta’. Selects the encoding for the time of flights.

multi_objective (bool): when True constructs a multiobjective problem (dv, T).

orbit_insertion (bool): when True the arrival dv is computed as that required to acquire a target orbit defined by e_target and rp_target.

e_target (float): if orbit_insertion is True this defines the target orbit eccentricity around the final planet.

rp_target (float): if orbit_insertion is True this defines the target orbit pericenter around the final planet (in m).

eta_bounds (list [float]): lower and upper bounds for the eta variable (defining the position of the DSM).

rp_ub (float): upper bound for the pericenter radius of the flybys (in planetary radii).

plot(x, ax=None, units=pk.AU, N=60, c_orbit='dimgray', c_lambert='indianred', c_ballistic='royalblue', leg_ids=[], figsize=(5, 5), **kwargs)

Plots the trajectory encoded into x in 3D axes.

Args:

x (list): The decision vector in the correct tof encoding.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): the 3D axis. If None a new one is created.

units (float, optional): length units to be used. Defaults to pk.AU.

N (int, optional): number of points to be used in the plots. Defaults to 60.

c_orbit (str, optional): color of the planetary orbits. Defaults to ‘dimgray’.

c (str, optional): color of the trajectory. Defaults to ‘indianred’.

figsize (tuple, optional): size of the figure. Defaults to (5, 5).

**kwargs: Additional keyword arguments to pass to the trajectory plotting functions.

pretty(x)

• x: encoded trajectory

Prints human readable information on the trajectory represented by the decision vector x

Example:

print(prob.pretty(x))

class pykep.trajopt.pl2pl_N_impulses(start='earth', target='venus', N_max=3, tof=[20., 400.], vinf=[0., 4.], phase_free=True, multi_objective=False, t0=None)

This class is a pygmo (http://esa.github.io/pygmo/) UDP representing a single leg transfer between two planets allowing up to a maximum number of impulsive Deep Space Maneuvers.

The decision vector is:

[t0,T] + [alpha,u,v,V_inf] * (N-2) + [alpha] + ([tf])

… in the units: [mjd2000, days] + [nd, nd, m/sec, nd] + [nd] + [mjd2000]

Each time-of-flight can be decoded as follows, T_n = T log(alpha_n) / sum_i(log(alpha_i))

Note

The resulting problem is box-bounded (unconstrained). The resulting trajectory is time-bounded.

Args:

start (planet): initial body.

target (planet): target body.

N_max (int): maximum number of impulses.

tof_bounds (list): the box bounds [lower,upper] for the total time of flight (days).

DV_max_bounds (list): the box bounds [lower,upper] for each DV magnitude (km/sec).

phase_free (bool): when True, no rendezvous condition are enforced and start and arrival anomalies will be free.

multi_objective (bool): when True, a multi-objective problem is constructed with DV and time of flight as objectives.

t0_bounds (list): the box bounds on the launch window containing two pykep.epoch. This is not needed if phase_free is True.

plot(x, ax=None, units=pk.AU, N=60, c_orbit='dimgray', c_segments=['royalblue', 'indianred'], figsize=(5, 5), **kwargs)

Plots the trajectory encoded into x in 3D axes.

Args:

x (list): The decision vector in the correct tof encoding.

ax (mpl_toolkits.mplot3d.axes3d.Axes3D, optional): The 3D axis to plot on. Defaults to None.

units (float, optional): The unit scale for the plot. Defaults to pk.AU.

N (int, optional): The number of points to use when plotting the trajectory. Defaults to 60.

c_orbit (str, optional): The color of the planet orbits. Defaults to ‘dimgray’.

c_segments (list, optional): The colors to alternate the various trajectory segments (inbetween DSMs). Defaults to [“royalblue”, “indianred”].

figsize (tuple): The figure size (only used if a*ax* is None and axis have to be created.), Defaults to (5, 5).

**kwargs: Additional keyword arguments to pass to the trajectory plot (common to Lambert arcs and ballistic arcs)

Returns:

mpl_toolkits.mplot3d.axes3d.Axes3D: The 3D axis where the trajectory was plotted.

plot_primer_vector(x, N=200, ax=None)

Plots the primer vector magnitude along the trajectory encoded in x.

Args:

x (list): The decision vector in the correct tof encoding.

N (int, optional): The number of points to use when plotting the primer vector. Defaults to 200.

ax (matplotlib.axes.Axes, optional): The axis to plot on. Defaults to None.

Returns:

matplotlib.axes.Axes: The axis where the primer vector was plotted. tuple: A tuple containing the grid and the primer vector magnitude.

Utilities

In order to facilitate the use of the classes in this module, some utilities are provided.

---

pykep.trajopt.primer_vector(DVi, DVj, Mji, Mjk)

This function computes the primer vector in a point k, relative to finite impulses in i and j.

Args:

DVi (ndarray - (3,)): the impulse at point i.

DVj (ndarray - (3,)): the impulse at point j.

Mji (ndarray - (6,6)): the state transition matrix from i to j (dxj = Mji dxi).

Mjk (ndarray - (6,6)): the state transition matrix from k to j (dxj = Mji dxi).

Returns:

tuple: The primer vector, the Aik matrix, the Ajk matrix.

Note:

The impulse transfer matrix Anm is defined as those matrices that allow to compute the variation of the impulse at point n given the variation of the impulse at point m. In formal terms, dDVn = Anm dDVm. All variations are such that the terminal state is kept fixed.

pykep.trajopt.primer_vector_surrogate(DVk, Mki, Mkj)

This function computes the surrogate primer vector at two points i and j, corresponding to a single finite impulse at point k.

Args:

DVk (ndarray - (3,)): the impulse at point k.

Mki (ndarray - (6,6)): the state transition matrix from i to k (dxk = Mki dxi).

Mkj (ndarray - (6,6)): the state transition matrix from j to k (dxk = Mkj dxj).

Returns:

tuple: The surrogate primer vector, the Aij matrix, the Akj matrix.

Note:

The impulse transfer matrix Anm is defined as that matrix that allows to compute (in this surrogate case) the variation of the impulse at point n given the variation of the impulse at point m. In formal terms, dDVn = Anm dDVm. All variations are such that the terminal state is kept fixed.

class pykep.trajopt._launchers

This class contains an API to access the tabulated data on the mass launchers can deliver to a certain declination / vinf.

Note

In pykep the object pykep.trajopt.launchers is an instance of this class created upon import and is used to access launchers data.

Examples:

>>> import pykep as pk
>>> mass = pk.trajopt.launchers.soyuzf(4.5, 33.21)
>>> mass2 = pk.trajopt.launchers.atlas501(4.5, 33.21)

ariane5(vinfs, decls)

Estimates the mass that the Ariane5 launcher can deliver to a certain vinf and declination, assuming a launch from Kourou. Data provided to ESOC by Arianespace when Ariane launch for ExoMars was an option. If the inputs are arrays, then a mesh is constructed and the mass is returned on all points of the mesh.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

decls (float or numpy.ndarray (N)):: the declinations in degrees (in the ECF frame).

Returns:

numpy.ndarray (N, N): masses delivered to escape with said declinations and magnitudes.

atlas501(vinfs, decls)

Estimates the mass that the Atlas 501 launcher can deliver to a certain vinf and declination. If the inputs are arrays, then a mesh is constructed and the mass is returned on all points of the mesh.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

decls (float or numpy.ndarray (N)):: the declinations in degrees (in the ECF frame).

Returns:

numpy.ndarray (N, N): masses delivered to escape with said declinations and magnitudes.

atlas551(vinfs)

Estimates the mass that the Atlas 551 launcher can deliver to a certain vinf

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

Returns:

numpy.ndarray (N): masses delivered to escape with said vinf magnitudes.

soyuzf(vinfs, decls)

Estimates the mass that the Soyutz-Fregat launcher can deliver to a certain vinf and declination. If the inputs are arrays, then a mesh is constructed and the mass is returned on all points of the mesh.

Args:

vinfs (float or numpy.ndarray (N)): the hyperbolic escape velocities in km/s.

decls (float or numpy.ndarray (N)):: the declinations in degrees (in the ECF frame).

Returns:

numpy.ndarray (N, N): masses delivered to escape with said declinations and magnitudes.