PyKEP: Open-Source Toolkit for Interplanetary Design

• PyKEP is an open-source software environment that integrates fast astrodynamics primitives such as Lambert solvers, Taylor propagators, and Sims–Flanagan low-thrust methods.
• It enables rapid assembly and evaluation of trajectory building blocks through methods like beam search, tree search, and porkchop grid generation for preliminary mission planning.
• The toolkit supports complex applications including multi-asteroid rendezvous, phasing value computations using orbital indicators, and efficient node ranking for global trajectory exploration.

Searching arXiv for the specified paper to ground the response. PyKEP is an open-source library focused on preliminary interplanetary mission design. In the treatment centered on "Designing Complex Interplanetary Trajectories for the Global Trajectory Optimization Competitions" (Izzo et al., 2015), it is presented as a software environment that wraps fast, well-tested astrodynamics primitives and transcriptions to quickly assemble trajectory building blocks, evaluate feasibility, and feed optimization or tree-search algorithms. The chapter uses PyKEP both to introduce fundamental problems of space flight mechanics and to formulate a multiple asteroid rendezvous strategy related to the 7th edition of the Global Trajectory Optimization Competition, with particular emphasis on Lambert transfers, low-thrust transcriptions, phasing indicators, clustering, and beam-search-based global exploration (Izzo et al., 2015).

1. Position within preliminary interplanetary design

PyKEP is used for a typical preliminary-design workflow in which bodies, ephemerides, and epochs are first defined; transfer legs are then computed via Lambert methods under impulsive assumptions or via transcribed low-thrust arcs; $\Delta v$, mass, and time budgets are sized; coarse grids such as porkchops are generated to identify windows; and a global search is performed using beam search, tree search, or heuristics guided by surrogate metrics, after which promising sequences are refined (Izzo et al., 2015).

Within that workflow, the chapter employs PyKEP to query ephemerides and simple planet models, including JPL low-precision models for Earth and Mars; solve the Kepler initial value problem with Lagrange coefficients through propagate_lagrangian; solve Lambert’s boundary value problem through lambert_problem, including multi-revolution branches when desired; numerically propagate constant-thrust arcs with a Taylor-series integrator through propagate_taylor; transcribe and evaluate low-thrust optimal control problems with the Sims–Flanagan method through sims_flanagan spacecraft and leg objects; compute and use phasing indicators and clustering; and orchestrate tree searches, including beam search and MOBS, to generate, rank, and prune trajectory candidates (Izzo et al., 2015).

The chapter also notes that PyKEP ships tools for gravity-assist sequence exploration, such as MGA and MGA-1DS, and that porkchop-style grid generation can be built by sweeping Lambert solves across departure and arrival grids. These are described as typical preliminary-design tools even though they are not the main focus of the chapter. A plausible implication is that the library is intended less as a monolithic mission optimizer than as a composable toolkit for rapid assembly of trajectory-analysis pipelines.

2. Core astrodynamics primitives and transcriptions

A central building block is Lambert’s problem, formulated as a two-body boundary-value problem,

$$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$

The inputs are $\mu$, $\mathbf r_s$, $\mathbf r_f$, $T=t_f-t_s$, the transfer geometry, and the multi-revolution index; the outputs are the departure and arrival velocities $v_s^\ast=v_1$ and $v_f^\ast=v_2$, with possibly multiple branches, including $0$-rev and higher-revolution solutions, and prograde or retrograde cases. PyKEP exposes this through the pattern $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$24 and uses Izzo’s solver. For zero-revolution transfers, the chapter reports near-constant runtime across geometries; in the chapter’s measurements, a $0$-rev Lambert problem took approximately $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$0 on a single thread of an Intel i7-4600U, corresponding to roughly $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$1k Lambert problems per second (Izzo et al., 2015).

For unperturbed two-body propagation, PyKEP wraps the Lagrange $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$2–$$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$3 formulation,

$$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$4

via propagate_lagrangian(rs, vs, dt, mu). The coefficients depend on $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$5, the universal variable $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$6, the Stumpff functions $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$7 and $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$8, and the geometry, allowing propagation without explicit numerical integration for the unperturbed two-body motion.

For constant-thrust propagation, the chapter uses propagate_taylor, described as a Taylor-series integrator for perturbed Keplerian motion with constant thrust and mass flow. For low-thrust optimal control problems, the principal transcription is the Sims–Flanagan method. The chapter’s representative pattern constructs a sims_flanagan.spacecraft, defines initial and final spacecraft states through sims_flanagan.sc_state, and builds a sims_flanagan.leg; mismatch constraints and throttle constraints are then evaluated through leg.mismatch_constraints() and leg.throttles_constraints(). This makes the transcription directly usable as a feasibility filter or as a component in broader search strategies (Izzo et al., 2015).

These components are used as interoperable primitives rather than as isolated solvers. Lambert arcs support impulsive screening and seeding, while Sims–Flanagan legs support low-thrust feasibility analysis and resource accounting. This suggests that PyKEP’s design emphasizes high-throughput evaluation of candidate transfers under preliminary-design assumptions.

3. Mechanical relations and trajectory evaluation

The chapter places PyKEP’s software abstractions within standard space-flight mechanics. In addition to Lambert’s boundary-value problem and the Kepler initial value problem,

$$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$9

it recalls elliptic Kepler’s equation,

$\mu$0

and the vis-viva relation,

$\mu$1

These are described as relations used internally by lambert_problem and ephemerides routines and as underpinning analytics in two-body propagation and orbital energy checks (Izzo et al., 2015).

Mass and velocity budgeting are organized around the Tsiolkovsky rocket equation,

$\mu$2

and, for constant thrust,

$\mu$3

For impulsive Lambert legs, $\mu$4 accounting is performed by matching encounter relative velocities:

$\mu$5

For low-thrust legs, mass flow enters the Sims–Flanagan and Taylor propagators directly.

The chapter also describes porkchop data generation as a grid over departure times $\mu$6 and arrival times $\mu$7, with metrics such as $\mu$8, $\mu$9, and time of flight. PyKEP facilitates these grids by looping ephemerides queries and Lambert solves. The concrete examples use planet.jpl_lp('earth'), planet.jpl_lp('mars'), and epochs in mjd2000, with explicit conversion of time of flight to seconds before calling lambert_problem. The text explicitly warns that misaligned units, for example mixing days and seconds, will corrupt time of flight and $\mathbf r_s$0; it also notes that the chapter’s examples often use non-dimensional units inside PyKEP and therefore require consistent unit handling when mixing m/s, AU/day, and related conventions (Izzo et al., 2015).

A common misconception in preliminary mission design is that coarse dynamical models are intrinsically inappropriate for algorithmic search. The chapter instead adopts heliocentric two-body dynamics, Keplerian asteroid ephemerides, and constant-thrust arcs for the preliminary phase. This does not assert sufficiency for final flight design; rather, it identifies a deliberate trade between dynamical fidelity and tractable global exploration.

4. GTOC7-inspired multiple-asteroid rendezvous formulation

The chapter’s principal case study is a multiple asteroid rendezvous problem related to GTOC7. The target set contains $\mathbf r_s$1 main-belt asteroids with Keplerian ephemerides. The spacecraft parameters are $\mathbf r_s$2, propellant mass $\mathbf r_s$3, $\mathbf r_s$4, and $\mathbf r_s$5. The mission duration is limited to $\mathbf r_s$6 years, with a minimum stay time per asteroid of $\mathbf r_s$7 days. The start date lies in $\mathbf r_s$8 MJD2000, and the spacecraft can be placed on a chosen starting asteroid at a chosen date. Rendezvous requires perfect position and velocity match; the dynamics are heliocentric two-body; and low-thrust propulsion is used for transfers, evaluated through Sims–Flanagan and Taylor propagation. The stated objective is to maximize the number of asteroid rendezvous under the remaining time and remaining mass resources (Izzo et al., 2015).

The global-search strategy is MOBS, a multi-objective beam search. A node is defined by the sequence of visited asteroids together with the non-dimensional remaining resources

$\mathbf r_s$9

At each node, the branch procedure first selects up to $\mathbf r_f$0 candidate next asteroids through $\mathbf r_f$1-nearest neighbors in the orbital-indicator metric $\mathbf r_f$2; for each candidate, it then solves an earliest-arrival low-thrust optimal control problem to obtain the earliest feasible arrival time $\mathbf r_f$3; if such a solution exists, it solves a set of fixed-arrival minimum-mass optimal control problems with

$\mathbf r_f$4

to capture feasible variants; and each feasible solution yields updated $\mathbf r_f$5 and $\mathbf r_f$6 and therefore an extended node (Izzo et al., 2015).

Nodes at the same depth are ranked through a node value $\mathbf r_f$7. The options tested in the chapter are best mass $\mathbf r_f$8, best time $\mathbf r_f$9, the soft minimum

$T=t_f-t_s$0

and Pareto dominance $T=t_f-t_s$1. The beam keeps the best $T=t_f-t_s$2 nodes for further expansion. The compact pseudocode in the chapter is:

$$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$25

The constraints handled in this formulation include departure and arrival windows $T=t_f-t_s$3 when computing phasing, mission-horizon tracking through $T=t_f-t_s$4 across legs, and enforcement of the minimum stay time $T=t_f-t_s$5 between rendezvous events. Maximum time of flight per leg is enforced implicitly by the earliest-arrival and fixed-arrival checks. The chapter also notes that launch $T=t_f-t_s$6 limits are not central to this belt-mission setup because the spacecraft is delivered to the starting asteroid (Izzo et al., 2015).

5. Phasing value, surrogates, and clustering

The chapter introduces the phasing value as a multi-objective measure of the transfer quality between asteroids. With the spacecraft $T=t_f-t_s$7 at asteroid $T=t_f-t_s$8 at time $T=t_f-t_s$9, Lambert problems are considered to a target asteroid $v_s^\ast=v_1$0 over $v_s^\ast=v_1$1, with thrust feasibility constrained by

$v_s^\ast=v_1$2

The associated objectives are

$v_s^\ast=v_1$3

subject to $v_s^\ast=v_1$4. The phasing value is the hypervolume of the Pareto front for $v_s^\ast=v_1$5, using the reference point

$v_s^\ast=v_1$6

A larger hypervolume corresponds to better phasing, meaning more and deeper feasible trade-offs between lower $v_s^\ast=v_1$7 and earlier arrival (Izzo et al., 2015).

Because direct phasing-value computation is expensive, the chapter introduces two fast surrogates. The Euclidean indicator is

$v_s^\ast=v_1$8

which captures proximity in phase space but cannot distinguish converging from diverging relative motion. The orbital indicator is based on a linearized transfer model:

$v_s^\ast=v_1$9

with

$v_f^\ast=v_2$0

If asteroid $v_f^\ast=v_2$1 at epoch $v_f^\ast=v_2$2 is mapped to

$v_f^\ast=v_2$3

then $v_f^\ast=v_2$4, i.e. Euclidean distance in $v_f^\ast=v_2$5D (Izzo et al., 2015).

This representation enables $v_f^\ast=v_2$6-d-tree acceleration. The chapter gives complexity figures of $v_f^\ast=v_2$7 to build the tree and $v_f^\ast=v_2$8 to query $v_f^\ast=v_2$9 nearest neighbors. In the reported tests, one $0$0-NN computation, including tree construction, took approximately $0$1 s, whereas computing all phasing values $0$2 took about $0$3 minutes; Kendall-$0$4 correlations confirmed that $0$5 performed better than $0$6 as a ranking surrogate (Izzo et al., 2015). The practical implication is direct: a computationally expensive multi-objective transfer-quality calculation can be replaced by an inexpensive metric that preserves the search order sufficiently well for large-scale branching.

The same orbital-indicator vectors support clustering. A cluster is defined as a subset of asteroids at a given epoch $0$7 whose pairwise $0$8 distances are small, specifically within $0$9. The chapter proposes DBSCAN with Euclidean metric on the orbital-indicator vectors, using parameters $0$0 and $0$1, where $0$2 is interpretable as a $0$3 budget under the linear model and $0$4 is the minimum number of neighbors for a core point. The stated procedure is to compute the $0$5D vectors, run DBSCAN, identify core and border points, discard outliers, and then choose the start epoch by scanning for epochs with large clusters, selecting the initial asteroid from the core of the largest cluster. Runtime is described as scaling well for $0$6–$0$7 because the index build is $0$8 and region queries are fast (Izzo et al., 2015).

6. Reported performance, limitations, and comparative context

Under the parameter set $0$9, $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$00, $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$01, $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$02, $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$03 years, and $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$04 days, the chapter reports that MOBS consistently found sequences of length $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$05, occasionally found $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$06, and did not reach $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$07 in the reported runs; the best sequence length reported in the competition was $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$08 (Izzo et al., 2015). For node ranking, best time $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$09 was the most effective node value, while best mass $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$10 performed poorly; soft-min and Pareto-based ranking were intermediate. The chapter interprets this as evidence that time is the scarcer resource in the stated setup.

Clustering materially affected success rates. Starting inside a large cluster with $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$11 and $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$12 significantly improved outcomes, and epochs near approximately $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$13 MJD2000 exhibited notably large clusters, identified as good-phasing epochs (Izzo et al., 2015). This suggests that start-state selection is not merely an initialization detail but a dominant determinant of reachable sequence length.

The computational campaign was substantial. The chapter reports approximately $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$14 million optimal control problem evaluations across $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$15 MOBS runs; one MOBS run required approximately $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$16 hour on a $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$17 GHz multi-core machine; and the overall campaign ran in about $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$18 days on $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$19 cores with $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$20 hyperthreads (Izzo et al., 2015). Several practical guidelines are stated explicitly: keep mr=0 for most Lambert calls unless multi-revolution branches are needed; use $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$21-d trees for nearest-neighbor selection and DBSCAN neighborhood operations; use the Taylor integrator for constant-thrust propagation because its runtime scales roughly linearly with integration time; enforce throttle constraints in Sims–Flanagan legs; and drive mismatch constraints near zero for feasibility.

The limitations are also explicit. Full phasing-value computation through a multi-objective Lambert search with hypervolume is expensive at scale, whereas the orbital indicator $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$22 provides an effective surrogate with strong rank correlation. Global search trees grow rapidly, so beam search with cluster-guided starts is used to provide predictable runtime per depth, approximately $$\ddot{\mathbf r} = -\mu \frac{\mathbf r}{r^3},\quad \mathbf r(t_s) = \mathbf r_s,\quad \mathbf r(t_f) = \mathbf r_f.$$23 cost per branch, and thereby budget control (Izzo et al., 2015).

In comparative terms, PyKEP is described as providing fast, open-source implementations of Lambert solvers, ephemerides, low-thrust transcriptions such as Sims–Flanagan, and numerical propagators, making rapid prototyping of global search strategies tractable and allowing integration with external optimizers including IPOPT, SNOPT, and WORHP without vendor lock-in. Tree-search strategies such as beam search and Monte Carlo Tree Search are characterized as algorithmic layers on top of these building blocks. Compared with dedicated flight-dynamics environments, PyKEP is said to trade high-fidelity perturbations for speed and simplicity in preliminary design, precisely the regime exploited in the chapter through Keplerian asteroids, a heliocentric two-body baseline, and constant-thrust arcs (Izzo et al., 2015).