The pontryagin module

This module contains the class pykep.pontryagin.leg, which allows one to efficiently construct low-thrust trajectories using the indirect optimal control transcription alla moda di Pontryagin’s maximum principle. pykep.pontryagin.leg transcribes the two-point boundary value problem resulting from Pontryagin’s maximum principle , in which one must correctly choose (via and optimiser, e.g. SNOPT) the nondimensional costate variables (nondimensional variables each corresponding to the respective components of the state) to lead the dynamical system (i.e. spacecraft dynamics) to chosen boundary conditions (within some tolerance). The user can choose the parameters type(freemass) == bool and type(freetime) == bool to enforce transversality conditions on the arrival mass and time. If freemass == True than the final mass may vary; if freetime == True than the final time may vary.

This configurability allows one to correctly construct any of the following:

  • single shooting trajectory optimisation problems

  • multiple shooting trajectory optimisation problem

  • trajectory optimisation problems with flybys

The list of classes and the detailed documentation follows:

Name

Type

Description

pykep.pontryagin.leg

class

represents one leg in the Pontryagin model

Detailed Documentation

class pykep.pontryagin.leg(object)

Indirect optimal control transcription trajectory leg.

This class represents an indirect optimal control transcription (alla moda di Pontryagin’s maximum principle) of a generic two-point boundary trajectory.

Attributes:

  • t0 (float): Departure time [mjd2000].

  • x0 (numpy.ndarray): Cartesian departure state [m, m, m, m/s, m/s, m/s, kg].

  • l0 (numpy.ndarray): Costate variables [ND, ND, ND, ND, ND, ND, ND].

  • tf (float): Arrival time [mjd2000].

  • xf (numpy.ndarray): Cartesian arrival state [m, m, m, m/s, m/s, m/s, kg].

  • sc (pykep.sims_flanagan.spacecraft): Generic spacecraft with propulsive properties.

  • mu (float): Gravitational parametre of primary body [m^3/s^2].

  • freemass (bool): Activates final mass transversality condition.

  • freetime (bool): Activates final time transversality condition. Allows final time to vary.

  • alpha (float): Homotopy parametre.

  • bound (bool): Activates bounded control.

  • nec (int): Number of equality constraints.

  • trajectory (numpy.ndarray): Array of nondimensional fullstates.

  • times (numpy.ndarray): Vector of nondimensional integration times.

__init__(t0=None, x0=None, l0=None, tf=None, xf=None, sc=pykep.sims_flanagan.spacecraft(1000, 0.3, 2500), mu=pykep.MU_SUN, freemass=True, freetime=True, alpha=1, bound=True)

Initialises an indirect optimal control transcription trajectory leg.

Args:

  • t0 (pykep.epoch, None): Departure time [mjd2000].

  • x0 (pykep.sims_flanagan.sc_state, None): Cartesian departure state [m, m, m, m/s, m/s, m/s, kg].

  • l0 (numpy.ndarray, list, tuple, None): Costate variables [ND, ND, ND, ND, ND, ND, ND].

  • tf (pykep.epoch, None): Arrival time [mjd2000].

  • xf (pykep.sims_flanagan.sc_state, None): Cartesian arrival state [m, m, m, m/s, m/s, m/s, kg].

  • sc (pykep.sims_flanagan.spacecraft): Generic spacecraft with propulsive properties.

  • mu (float, int): Gravitational parametre of primary body [m^3/s^2].

  • freemass (bool): Activates final mass transversality condition.

  • freetime (bool): Activates final time transversality condition. Allows final time to vary.

  • alpha (float, int): Homotopy parametre, governing the degree to which the theoretical control law is intended to reduce propellant expenditure or energy. Setting the parametre to 1 enforces a mass-optimal control law, with a characteristic bang-bang control profile (either full throttle or off). Setting the parametre to 0 enforces a pure quadratic control law.

  • bound (bool): Activates bounded control, in which the control throttle is bounded between 0 and 1, otherwise the control throttle is allowed to unbounded.

Raises:

  • TypeError: If sc is not supplied as an instance of pykep.sims_flanagan.spacecraft.

  • TypeError: If mu is neither supplied as an instance of float nor int.

  • ValueError: If mu is not supplied as a positive number.

  • TypeError: If either freemass, freetime, or bound is not supplied as an instance of bool.

  • AttributeError: If equality constraint dimension cannot be determined form freemass and freetime.

  • TypeError: If alpha is neither supplied as an instance of float nor int.

  • ValueError: If alpha is not in between 0 and 1.

  • ValueError: If alpha == 1 and bound == False. Control cannot be unbounded with mass-optimal control.

  • TypeError: If either t0 or tf is not supplied as an instance of pykep.epoch.

  • ValueError: If departure time, t0.mjd2000 does not occur before arrival time, t`f.mjd2000`.

  • TypeError: If either x0 or xf is not supplied as an instance of pykep.sims_flanagan.sc_state.

  • TypeError: If costate, l0, is neither supplied as an instance of numpy.ndarray, list, nor tuple.

  • TypeError: If costate, l0, is does not have 7 elements. Each element of l0 corresponds to the respective elements of x0.

Examples:
>>> l = pykep.pontryagin.leg()
>>> l = pykep.pontryagin.leg(t0, x0, l0, tf, xf)

Note

Typically spacecraft trajectories are optimised to reduce final mass, which theoretically results in a bang-bang control profile. However, a bang-bang control profile (either 0 or 1) is problematic to an optimiser due to its discontinous nature. The trajectory optimisation problem becomes easier when one first optimises with unbounded quadratic control (alpha == 0 and bound == False), stores the solution, then uses the solution to optimise with bounded quadratic control (alpha == 0 and bound == True). Then using the solution from bounded quadratic control, one can incrementally optimise for increasing homotopy parametres (e.g. alphas = numpy.linspace(0, 1, 200)) until a mass-optimal control solution converges.

Note

Boundary conditions t0, x0, l0, tf, and xf are optional in the constructor. If some, but not all, boundary conditions are supplied, they are not set and disregarded. If the boundary conditions are not set upon construction (__init__), they must be set with set(t0, x0, l0, tf, xf).

Note

mismatch_constraints will not work unless the boundary conditions have been set either through __init__ or set.

set(t0, x0, l0, tf, xf)

Sets the departure and arrival boundary conditions of the trajectory.

This instance method sets the departure time t0, departure state x0, departure costate l0, arrival time tf, arrival state xf.

Args:

  • t0 (pykep.epoch): Departure time [mjd2000].

  • x0 (pykep.sims_flanagan.sc_state): Cartesian departure state [m, m, m, m/s, m/s, m/s, kg].

  • l0 (numpy.ndarray, list, tuple): Costate variables [ND, ND, ND, ND, ND, ND, ND].

  • tf (pykep.epoch): Arrival time [mjd2000].

  • xf (pykep.epoch): Cartesian arrival state [m, m, m, m/s, m/s, m/s, kg].

Raises:

  • TypeError: If either t0 or tf is not supplied as an instance of pykep.epoch.

  • TypeError: If either x0 or xf is not supllied as an instance of pykep.sims_flanagan.sc_state.

  • TypeError: If l0 is neither supllied as a numpy.ndarray, list, nor tuple.

Note

If args were not supplied in __init__, set must be called before calling mismatch_constraints.

Examples:
>>> l = leg() # t0, x0, l0, tf, xf not set yet
>>> l.set(t0, x0, l0, tf, xf) # now they're set
mismatch_constraints(atol=1e-05, rtol=1e-05)

Returns the nondimensional mismatch equality constraints of the arrival boundary conditions.

This method propagates the nondimensional dynamics of the spacecraft from the departure time t0 to the arrival time tf, then evaluates the nondimensional mismatch between the desired arrival state xf, arrival mass costate lmf == 0 (if freemass == True), and arrival Hamiltonian H == 0 (if freetime == True) and the supplied desired arrival boundary conditions.

Args:

  • atol (float, int): Absolute integration solution tolerance.

  • rtol (float, int): Relative integration solution tolerance.

Returns:

  • numpy.ndarray: Equality constraints vector ceq composed of the arrival mismatch in position & velocity, arrival mass costate if freemass == True, arrival mismatch in mass if freemass == False, and arrival Hamiltonian if freetime == True.

ceq = [drf, dvf, dmf]    # freemass == False; freetime == False
ceq = [drf, dvf, lmf]    # freemass == True; freetime == False
ceq = [drf, dvf, lmf, H] # freemass == True; freetime == True
Raises:

  • TypeError: If either atol or rtol is supplied as neither an instance of float nor int.

  • AttributeError: If boundary conditions t0, x0, l0, tf, and x0 have not been set through either __init__ or set.

Note

This method uses an explicit Runge-Kutta method of order (4)5 due to Dormand & Prince with adaptive stepsize control. Smaller values of atol and rtol will increase the accuracy of a converged trajectory optimisation solution. However, smaller values result in slower integration executions and thus slower optimiser iterations. It is preferable to first solve a trajectory optimisation problem with higher values of atol and rtol to quickly converge to an approximate solution, subsequently resolving the problem with smaller tolerance values using the previous solution. This method is akin to mesh refinement in direct trajectory optimisation.

Note

State arrival boundary conditions in position and velocity are always returned. If the final mass transversality condition is activated (freemass == True), than the final mass costate is returned as part of the equality constraint vector. Otherwise (freemass == False), the final mass mismatch is returned as part of the equality constraint vector. If the final time transversality condition is activated freetime == True, then than the final Hamiltonian will be returned as part of the equality constraint vector. Otherwise (freetime == False), the Hamiltonian will not be supplied as part of the equality constraint vector.

Examples:
>>> l = pykep.pontryagin.leg(freemass=True, freetime=True)
>>> l.set(t0, x0, l0, tf, xf)
>>> l.mismatch_constraints(atol=1e-10, rtol=1e-10)
[ -1.15402379e+01  -2.23345886e+00  -9.30022917e-01  -2.53440778e+00
  -3.44246359e+00  -3.96669697e-01  -2.82967399e+03  -8.58526037e-01]
>>> l.nec
8

>>> l = pykep.pontryagin.leg(freemass=False, freetime=True)
>>> l.set(t0, x0, l0, tf, xf)
>>> l.mismatch_constraints(atol=1e-10, rtol=1e-10)
[-0.82435194  0.4375439   0.04615264  0.69192818 -0.18327442 -0.00930848
  0.46335841  1.66752446]
>>> l.nec
8

>>> l = pykep.pontryagin.leg(freemass=True, freetime=False)
>>> l.set(t0, x0, l0, tf, xf)
>>> l.mismatch_constraints(atol=1e-10, rtol=1e-10)
[  4.95782514e+00  -7.94403974e+00  -1.10158930e-02   7.27278513e+00
  -1.70998019e+00   9.13925064e-02  -1.45194673e+06]
>>> l.nec
7

>>> l = pykep.pontryagin.leg(freemass=False, freetime=False)
>>> l.set(t0, x0, l0, tf, xf)
>>> l.mismatch_constraints(atol=1e-10, rtol=1e-10)
[-0.90868565  0.23238016  0.04596579  0.61543688 -0.50023407 -0.0058185
  0.62170522]
>>> l.nec
7

>>> l = pykep.pontryagin.leg()
>>> l.mismatch_constraints(atol=1e-10, rtol=1e-10)
AttributeError: Cannot propagate dynamics, as boundary conditions t0, x0, l0, tf, and xf have not been set. Use set(t0, x0, l0, tf, xf) to set boundary conditions.
get_states(atol=1e-12, rtol=1e-12)

Returns the trajectory data of time, states, costates, and controls.

This method propagates the spacecrafts dynamics with a chosen atol rtol from t0, x0, l0 to t0, then returns a numpy.ndarray with the number of rows corresponding to the number of data points along the trajectory (as determinded by atol and rtol), and 19 columns.

Args:

  • atol (float, int): Absolute integration solution tolerance.

  • rtol (float, int): Relative integration solution tolerance.

Returns:

  • numpy.ndarray: Trajectory data array of shape (npts, 20), where npts is the number of data points along the trajectory, as determind by the integration error tolerances atol and rtol.

The returned array is characterised by:

[[t0, x0, y0, z0, vx0, vy0, vz0, m0, lx0, ly0, lz0, lvx0, lvy0, lvz0, lm0, obj, u0, ux0, uy0, uz0, H0],
...
 [tf, xf, yf, zf, vxf, vyf, vzf, mf, lxf, lyf, lzf, lvxf, lvyf, lvzf, lmf, obj, uf, uxf, uyf, uzf, Hf]]
Raises:

  • TypeError: If either atol or rtol is supplied as neither an instance of float nor int.

  • AttributeError: If boundary conditions t0, x0, l0, tf, and x0 have not been set through either __init__ or set.

Note

The returned array has the units of [mjd2000, m, m, m, m/s, m/s, m/s, kg, ND, ND, ND, ND, ND, ND, ND, s, ND, ND, ND, ND, ND].

Note

Setting either atol or rtol smaller than 1e-12 may result in numerical integration difficulties.

Examples:
>>> l = leg(t0, x0, l0, tf, xf)
>>> traj = l.get_states(atol=1e-12, rtol=1e-12)
>>> t = traj[:, 0] # times
>>> r = traj[:, 1:4] # positions
>>> v = traj[:, 4:7] # velocities
>>> m = traj[:, 7] # masses
>>> lm = traj[:, 14] # mass costates
>>> u = traj[:, 15] # control throttles
>>> H = traj[:, 19] # Hamiltonians
>>> plt.plot(t, u) # plot the control throttle history
>>> plt.show()
plot_traj(axis, mark='k.-', atol=1e-11, rtol=1e-11, units=pykep.AU)

Plots trajectory onto a 3D axis.

Args:

  • axes (matplotlib.axes._subplots.Axes3DSubplot): 3D axis onto which to plot the trajectory.

  • mark (str): Marker style.

  • atol (float, int): Absolute integration solution tolerance.

  • rtol (float, int): Relative integration solution tolerance.

  • units (float, int): Length unit by which to normalise data.

  • quiver (bool): Activates the visualization of the throttle arrows

  • length (float): Length of thrust arrow if quiver is True

Raises:

  • TypeError: If axes is not an instance of mpl_toolkits.mplot3d.Axes3D.

  • TypeError: If mark is not an instance of str.

  • TypeError: If either atol, rtol, or units is neither an instance of float nor int.

Examples:
>>> sc = pk.sims_flanagan.spacecraft(1000, 0.3, 2500) # spacecraft
>>> p0 = pk.planet.jpl_lp("earth")
>>> pf = pk.planet.jpl_lp("mars")
>>> t0 = pk.epoch(0)
>>> tf = pk.epoch(1000)
>>> r0, v0 = p0.eph(t0)
>>> rf, vf = pf.eph(tf)
>>> x0 = pk.sims_flanagan.sc_state(r0, v0, sc.mass)
>>> xf = pk.sims_flanagan.sc_state(rf, vf, sc.mass/10)
>>> l0 = np.random.randn(7)
>>> l = leg(t0, x0, l0, tf, xf)
>>> fig = plt.figure()
>>> axes = fig.gca(projection='3d')
>>> l.plot(axes)
>>> plt.show()
plot(x, y, mark='k.-', atol=1e-12, rtol=1e-12, unitsx=1, unitsy=1, xlabel=False, ylabel=False)

Plots in two dimensions of the leg’s trajectory data.

keys = ['t', 'x', 'y', 'z', 'vx', 'vy', 'vz', 'm', 'lx', 'ly', 'lz', 'lvx', 'lvy', 'lvz', 'lm', 'u', 'ux', 'uy', 'uz', 'H']
Args:

  • x (str): x-axis dimension in keys.

  • y (str): y-axis dimension in keys.

  • mark (str): Marker style.

  • atol (float, int): Absolute integration solution tolerance.

  • rtol (float, int): Relative integration solution tolerance.

  • unitsx (float, int): Unit by which to normalise x-axis data.

  • unitsy (float, int): Unit by which to normalise y-axis data.

  • xlabel (str, bool): x-axis label. If label is False, no label is placed; if True, dimension name is placed.

  • ylabel (str, bool): y-axis label. If label is False, no label is placed; if True, dimension name is placed.