Ascent Trajectory Optimizer

Gravity turn ascent profile optimization with a staged delta-v budget model.

A launch vehicle converts propellant into orbital speed, and most of what it produces never reaches the orbit. This library flies a staged vehicle from the pad to a circular orbit over a rotating spherical Earth, reports where every metre per second went, and then searches the guidance parameters for the profile that carries the most payload without breaking the dynamic pressure or axial load limits. The delta-v budget is the deliverable; the optimiser exists to produce one worth reading.

Results

Every number below is printed by one of the commands in Running it, or is a ratio of two numbers that are. The configuration is the default one throughout: the Falcon 9 Block 5 in its expendable configuration, a 200 km circular target orbit at 28.5 degrees inclination, a 0.5 s integration step, a 35 kPa dynamic pressure limit and a 6 g axial load factor limit.

The budget, and the proof that it closes

From examples/02_optimize_payload.py.

Term	m/s
Gravity loss	1528.3
Drag loss	48.4
Steering loss	66.3
Earth rotation credit (centrifugal term)	-4.2
Total losses	1638.8
Delta-v delivered (Earth-relative speed gain)	7362.7
Ideal delta-v (integral of thrust over mass)	9001.5
Closure residual	9.32e-07

The propulsion system produces 9001.5 m/s. The vehicle keeps 7362.7 m/s of it, which is the Earth-relative speed a 200 km circular orbit needs once the 421.6 m/s the rotating Earth contributes at that radius is taken off the 7784.3 m/s inertial circular speed. The other 1638.8 m/s is spent, 93 percent of it on gravity.

The last row is the reason to trust the rest. The five budget terms are not reconstructed from the trajectory afterwards; they are integrated alongside it, one accumulator per term, so the identity

ideal = delivered + gravity loss + drag loss + steering loss - rotation credit

is a live check on the whole right-hand side rather than on a post-processing routine. It closes to 9.32e-07 m/s in 9001.5 m/s, one part in 1e10, which is the truncation error of the integrator and not a modelling gap. tests/test_budget.py asserts that the residual falls at the order of the scheme when the step is halved, so the number stays a measurement rather than a coincidence.

Drag costs less than one percent of the budget, which is not because the atmosphere is cheap but because the dynamic pressure constraint keeps the vehicle out of the dense part of it. The sensitivity study below puts a number on that distinction.

The profile that produces it

Quantity	Value
Payload to target orbit	25713 kg
Target orbit	200 km circular
Pitch kick time	11.83 s
Pitch kick angle	1.730 deg
First stage cut-off time	152.3 s
Altitude at cut-off	69.2 km
Relative speed at cut-off	2617 m/s
Second stage burn time	374.0 s
Time of insertion	530.3 s
Peak dynamic pressure	34.99 kPa
Dynamic pressure limit	35.00 kPa
Peak axial load factor	5.163 g
Axial load factor limit	6.000 g
Insertion altitude error	-4.587e-02 m
Insertion radial velocity	-1.556e-04 m/s
Insertion speed error	-3.037e-05 m/s
Largest constraint violation	4.587e-05
Optimiser	Nelder-Mead then SLSQP
Converged	yes
Optimiser status	Optimization terminated successfully
SLSQP iterations	53
Trajectory evaluations	899

The solve converges, and insertion is met to 4.6 cm of altitude and 1.6e-04 m/s of radial velocity. The dynamic pressure constraint is active at the optimum and the load factor constraint is slack at 5.163 g, which is the expected shape of the trade: a shallower ascent spends less delta-v on gravity but stays low and fast for longer, so the best profile is the shallowest one the structural limit allows. That is what the figure at the top of this page shows, and it is the one thing in these results a table cannot: the trajectory does not merely respect the limit, it rides it.

Against the published capability

Quantity	Value
Published payload	22800 kg
Modelled payload	25713 kg
Payload discrepancy	+2913 kg
Payload discrepancy, relative	+12.8 %
Published liftoff mass	549054 kg
Modelled liftoff mass	572013 kg
Liftoff mass discrepancy, relative	+4.2 %
Published first stage sea level thrust	7.607 MN
Sea level thrust from the impulses	7.460 MN
Sea level thrust discrepancy, relative	-1.9 %
Published orbit	low Earth orbit, 28.5 degree inclination, expendable
Source	SpaceX, Falcon User's Guide, September 2021

The model over-predicts the payload by 12.8 percent, and it should. Four things a real vehicle spends are absent from it, all of them deliberately rather than by oversight:

• no flight performance reserve, the propellant a real mission holds back against dispersions and never plans to burn,
• no residual or unusable propellant, no pressurant and no attitude control propellant,
• no throttling, so the model flies a higher average thrust than the real vehicle does through max q and before cut-off,
• no out of plane steering, because the model is planar.

Each of those would reduce the payload. They are omitted rather than estimated, because an estimate chosen to close the gap would be tuning the model to the number it is supposed to be predicting.

The two rows that are not results are there for the same reason. The stage masses are third-party estimates rather than published data, and summing them with the modelled payload gives a liftoff mass 4.2 percent above the figure SpaceX publishes for this vehicle; since the payload elasticity to second stage dry mass is -0.155, an input wrong by a few percent moves the answer by a few hundred kilogrammes on its own. And the four published propulsion figures are not exactly self-consistent: taking the mass flow rate from the vacuum thrust and vacuum specific impulse, and the sea level thrust from the sea level specific impulse, gives 7.460 MN against a published 7.607 MN. The model uses the derived value so that all four published numbers are used consistently through one nozzle relation, and reports the 1.9 percent disagreement rather than hiding it.

What the atmosphere costs is the constraint, not the drag

From examples/04_sensitivity.py, which perturbs each parameter by 2 percent in each direction and re-optimises the whole ascent for every case. Elasticity is d log(payload) / d log(parameter), by central difference. All twelve perturbed solves converged.

Parameter	Low (kg)	High (kg)	Range (kg)	Elasticity
first stage vacuum specific impulse (2%)	24819	26485	+1666	+1.624
second stage vacuum specific impulse (2%)	24803	26611	+1808	+1.758
first stage dry mass (2%)	25783	25647	-136	-0.132
second stage dry mass (2%)	25793	25633	-160	-0.155
drag coefficient (2%)	25718	25707	-11	-0.010
dynamic pressure limit (2%)	25475	25740	+264	+0.258

This design is propulsion limited, not aerodynamics limited. Specific impulse dominates everything else by an order of magnitude, and the second stage value slightly beats the first because it acts on the whole remaining mass all the way to insertion. The drag coefficient is worth almost nothing: a 20 percent error in the drag data, far larger than the tabulated curve is likely to be wrong by, would move the payload by about 50 kg.

The dynamic pressure limit is worth roughly 24 times the drag coefficient. That is the same point stated the other way round, and it is the useful one: what the atmosphere costs this vehicle is not the drag it produces but the constraint it imposes on how early the vehicle may turn over. The response is also one-sided in a way the table hides. Tightening the limit by 2 percent costs 238 kg against the 25713 kg baseline, while loosening it by 2 percent buys only 27 kg, because a little way above the limit the constraint stops being active and the unconstrained optimum takes over.

What a coast is worth, and why the model now has one

From examples/05_high_orbit_coast.py, solving the same 600 km circular orbit twice, once with the second stage burning continuously to insertion and once with a burn, coast, burn sequence. Both runs converge from the default initial guess.

Quantity	Direct ascent	Burn, coast, burn
Target orbit	600 km circular	600 km circular
Payload to target orbit	16470 kg	24005 kg
Second stage first burn	374.0 s	369.8 s
Ballistic coast	none	3110.5 s
Second stage second burn	none	4.2 s
Time of insertion	530.3 s	3640.8 s
Peak dynamic pressure	32.22 kPa	35.00 kPa
Largest constraint violation	1.744e-05	6.711e-05
Converged	yes	yes
Payload gained by coasting		+7536 kg

A direct ascent to 600 km has to thrust its way up, and pays gravity loss for every second of the climb. Given a coast the optimiser finds the textbook alternative on its own, without being told what it is: burn onto a transfer ellipse, coast most of an hour to the top of it for free, and circularise with a 4.2 s burn. That is worth 7536 kg, or 46 percent of the direct ascent payload.

This capability closes a limitation that docs/design-notes.md previously recorded as open. The notes now record what closing it took, what it bought and what it cost, which is one more phase, two more decision variables and one more constraint. The 200 km result above is untouched by it: with a zero coast the phase schedule is the direct ascent exactly, and a test asserts so.

How the answer is produced

The ascent is integrated in Earth-centred inertial planar polar coordinates over a spherical, rotating Earth. Earth rotation enters the plane as an effective rate omega cos(i), exact for a due east launch from a site at latitude equal to the target inclination. The atmosphere is the U.S. Standard Atmosphere 1976, implemented below 86 km from its published lapse rates and layer base conditions rather than from its tabulated output, so the published table is an independent check on the implementation instead of its input. Drag is zero-lift against a Mach dependent drag coefficient. Thrust follows F(p) = F_vac - A_e p at constant mass flow, so the two published specific impulses of a stage fix its effective exit area.

Guidance is a gravity turn: a vertical rise, a pitch kick of a few degrees held over a short window, then a zero-lift turn with the thrust along the relative velocity vector, followed after staging by Lawden's linear tangent steering law. The linear tangent law is the exact optimal steering programme for the flat Earth, uniform gravity, vacuum, constant thrust problem, which makes it the natural two-parameter family for the vacuum phase rather than an arbitrary curve fit, and it is why the same law spans both burns when there is a coast between them.

The solve runs in two stages: a derivative free Nelder-Mead pass on a penalised objective, then SLSQP on the constrained problem. SLSQP suits the problem exactly, being a smooth objective with three nonlinear equality constraints for insertion, two nonlinear inequality constraints for the path limits and simple bounds in six variables, and its gradients are forward differences, which costs one trajectory integration per variable per iteration.

The first pass is not a precaution. A gravity turn is violently sensitive to the pitch kick, and a gradient method handed a raw initial guess takes its first finite difference at an impact point, where the terminal state says nothing about the optimum.

Two degrees of kick angle separate a trajectory that overshoots to 519 km from one that hits the ground. That is why the search box on the kick angle is narrow and why the solve starts without gradients.

The library is five layers, each with one job. Nothing in model or algorithm performs input or output, nothing outside analysis formats text or draws, and nothing in examples contains logic.

Layer	What lives there
model	Constants, the U.S. Standard Atmosphere, aerodynamics, gravity, the vehicle, and the Falcon 9 built from published specifications
algorithm	Guidance laws and optimiser wrappers, each behind a Protocol
pipeline	The phase schedule, the equations of motion with the budget accumulators, the fixed step integrator, and the maximum payload problem
analysis	The delta-v budget, the sensitivity study, the text tables and the figures
examples	Six scripts that build a problem, call the library, and print what comes back

The alternatives that were considered and rejected, including adaptive step integration, a global optimiser, trust region methods and an indirect formulation, are recorded with what each would have bought and cost in docs/design-notes.md.

Installation

Requires Python 3.12 or later.

git clone https://github.com/Eelis03/ascent-trajectory-optimizer.git
cd ascent-trajectory-optimizer
uv sync --all-extras --dev

Using pip instead of uv:

python -m venv .venv
.venv/bin/activate      # Windows: .venv\Scripts\activate
pip install -e ".[dev]"

Running it

import math

from ascent_optimizer import (
    FALCON_9_BLOCK_5,
    AscentProblem,
    TargetOrbit,
    delta_v_budget,
    format_delta_v_budget,
    solve_max_payload,
)

problem = AscentProblem(
    vehicle=FALCON_9_BLOCK_5,
    target=TargetOrbit(altitude_m=200_000.0, inclination_rad=math.radians(28.5)),
    time_step_s=0.5,
)
solution = solve_max_payload(problem)

print(f"converged: {solution.converged}")
print(f"payload: {solution.payload_mass_kg:.0f} kg")
print(format_delta_v_budget(delta_v_budget(solution.trace)))

Pass allow_coast=True to AscentProblem for the burn, coast, burn ascent. The decision vector grows from six variables to eight and the solve takes minutes rather than seconds, which is why it is not the default.

The example scripts produce every number on this page:

uv run python examples/01_reference_trajectory.py  # vehicle summary and kick angle sweep
uv run python examples/02_optimize_payload.py      # the budget and the published comparison
uv run python examples/03_step_convergence.py      # the order of the integrator
uv run python examples/04_sensitivity.py           # elasticities, twelve re-optimised solves
uv run python examples/05_high_orbit_coast.py      # direct against burn, coast, burn
uv run python examples/06_publication_figures.py   # the three figures in docs/figures

Every script that draws accepts --figure-dir, and every script that optimises accepts --max-iterations and --coarse-iterations to shorten the search. The two that solve for a high orbit or twelve perturbed cases take minutes; the rest take seconds.

The figures in docs/figures are snapshots of a single command,

uv run python examples/06_publication_figures.py

which re-solves the reference problem and overwrites all three. They are tracked in the repository rather than built on demand so that this page renders for a reader who never runs anything. Continuous integration runs the code that draws them, and the test suite checks that they exist, are real PNGs and stay inside the 250 kB budget, but it does not compare them byte for byte against the tracked copies: matplotlib output is not byte reproducible across platforms, versions or font configurations, and a check that fails for those reasons teaches a reader nothing.

How far the numbers can be trusted

The integrator is classical fourth order Runge-Kutta at a fixed step, with every phase divided into a whole number of equal substeps so that staging, jettison and each guidance switch land exactly on a step boundary. From examples/03_step_convergence.py, flying one fixed profile at halving steps:

Step	Final radius	Final relative speed	Budget closure residual
4.000 s	6658475.4688 m	7836.483858 m/s	1.828e-03 m/s
2.000 s	6663372.1680 m	7821.122123 m/s	1.355e-04 m/s
1.000 s	6663422.5762 m	7820.950868 m/s	1.168e-05 m/s
0.500 s	6663412.6613 m	7820.979552 m/s	8.369e-07 m/s
0.250 s	6663411.4905 m	7820.982761 m/s	5.556e-08 m/s
0.125 s	6663411.2328 m	7820.983372 m/s	3.567e-09 m/s

The closure residual falls by a factor of between 12 and 16 per halving, approaching 16 as the step falls, which is the fourth order behaviour of the scheme. The terminal state converges more slowly because the tabulated drag coefficient and the standard atmosphere lapse rate are continuous but not continuously differentiable, and contribute a lower order term wherever the trajectory crosses a table node. At the 0.5 s step used throughout, halving the step again moves the final radius by 1.17 m in 6663411 m, two orders inside the 1e-05 relative integration tolerance this project works to.

uv run pytest -q
uv run ruff check .
uv run mypy
uv run pytest --cov=src/ascent_optimizer --cov-report=term-missing

The suite is 136 tests in about 50 seconds and covers 98.75 percent of the package, 1123 statements with 14 missed. Coverage costs almost nothing here because the measurement runs on sys.monitoring rather than a per-line trace callback, which on code shaped like an integrator inner loop is the difference between 49 seconds and 182 seconds for the same 1123 statements and the same 14 missed. Continuous integration runs the suite on Linux and Windows with --cov-fail-under=96, that measurement rounded down and reduced by two, so a drop in coverage fails the build rather than being noticed later or not at all.

Three tiers: property and invariant tests over the mathematics, regression tests pinning recorded behaviour, and integration tests running every example script under a reduced iteration count. With gravity, drag, atmosphere and Earth rotation switched off and thrust held along the velocity vector, the integrator reproduces the Tsiolkovsky rocket equation to within 1e-06 m/s, which is 1.1e-10 of the value compared. The atmosphere reproduces all eight published layer base pressures to one unit in the last published digit.

Two rules govern every tolerance, and both are written into the docstrings of the tests they apply to.

Only reproducible values are pinned. A value produced by an iterative solve that did not converge is not reproducible, because the path an SLSQP run takes through a flat region near an optimum depends on the order in which a BLAS kernel reduces a dot product. Fixed-parameter trajectory outputs, which involve no search, are pinned at 1e-08 relative. The optimiser is pinned only as a converged objective, at 5 kg of payload, and convergence is asserted before any comparison is made. The parameter vector is never pinned.

Tolerances come from the measurement, not from the error that happened to be observed. The rocket equation tolerance is derived by fitting the fourth order error constant across three step sizes. The 1e-08 regression tolerance comes from counting the transcendental function evaluations along the trajectory, bounding the accumulated last-bit spread, and measuring how much the dynamics amplify an early perturbation. The 5 kg payload tolerance is twice the 2.5 kg spread measured across twelve converged runs at different ftol and finite difference settings, and remains an order of magnitude below the 48 kg that a change of one part in a thousand in second stage specific impulse produces, so it still fails on a real model change.

What is not modelled

The model is planar, gravity is a point mass with no J2 term, there is no throttling and no propellant reserve, no structural or thermal limit beyond dynamic pressure and axial load is computed, and the second stage steering is a two-parameter family rather than an optimal control, so the payload reported here is a lower bound on what an unrestricted control could reach. Several of the vehicle inputs are third-party estimates. Each of these is recorded in docs/design-notes.md with what removing it would take, alongside the one limitation that has been removed and what that cost.

A production ascent optimisation would not use a parameterised guidance family at all. It would use an indirect method or a pseudospectral direct method, and the design notes state what each buys and what each costs rather than merely naming them.

References

Models and methods:

• NOAA, NASA and USAF. "U.S. Standard Atmosphere, 1976". NOAA-S/T 76-1562, U.S. Government Printing Office, Washington DC, 1976. https://ntrs.nasa.gov/citations/19770009539
• National Imagery and Mapping Agency. "Department of Defense World Geodetic System 1984". NIMA TR8350.2, 3rd edition, amendment 1, 2000. Earth radius, gravitational parameter and rotation rate. https://earth-info.nga.mil/php/download.php?file=coord-wgs84
• Culler, G. J. and Fried, B. D. "Universal Gravity Turn Trajectories". Journal of Applied Physics, volume 28, number 6, 1957, pages 672 to 676. DOI 10.1063/1.1722828
• Lawden, D. F. "Optimal Trajectories for Space Navigation". Butterworths, London, 1963, chapter 3. The linear tangent steering law and the primer vector, which is why one steering programme spans a coast. https://catalog.hathitrust.org/Record/001620792
• McHenry, R. L., Long, A. D., Cockrell, B. F., Thibodeau, J. R. and Brand, T. J. "Space Shuttle Ascent Guidance, Navigation and Control". Journal of the Astronautical Sciences, volume 27, number 1, 1979, pages 1 to 38. https://ntrs.nasa.gov/citations/19790048206
• Sutton, G. P. and Biblarz, O. "Rocket Propulsion Elements". 9th edition, Wiley, 2017. ISBN 978-1-118-75388-0. Nozzle thrust against ambient pressure in section 3.3, launch vehicle drag coefficient and dynamic pressure in section 4.4.
• Hoerner, S. F. "Fluid-Dynamic Drag". Published by the author, Bakersfield CA, 1965, chapters 6 and 16. Drag of bodies of revolution, transonic and supersonic wave drag. https://catalog.hathitrust.org/Record/001621843
• Wiesel, W. E. "Spaceflight Dynamics". 3rd edition, Aphelion Press, 2010, chapter 4. The loss decomposition of a powered ascent. ISBN 978-1-4524-8938-1.
• Kraft, D. "A Software Package for Sequential Quadratic Programming". Technical Report DFVLR-FB 88-28, DLR German Aerospace Center, Institute for Flight Mechanics, Cologne, 1988. The SLSQP implementation SciPy wraps. https://degenerateconic.com/uploads/2018/03/DFVLR_FB_88_28.pdf
• Nelder, J. A. and Mead, R. "A Simplex Method for Function Minimization". The Computer Journal, volume 7, number 4, 1965, pages 308 to 313. DOI 10.1093/comjnl/7.4.308
• Hairer, E., Norsett, S. P. and Wanner, G. "Solving Ordinary Differential Equations I: Nonstiff Problems". 2nd edition, Springer, 1993, section II.1. DOI 10.1007/978-3-540-78862-1
• Betts, J. T. "Survey of Numerical Methods for Trajectory Optimization". Journal of Guidance, Control, and Dynamics, volume 21, number 2, 1998, pages 193 to 207. DOI 10.2514/2.4231
• Ross, I. M. and Karpenko, M. "A review of pseudospectral optimal control: from theory to flight". Annual Reviews in Control, volume 36, number 2, 2012, pages 182 to 197. DOI 10.1016/j.arcontrol.2012.09.002
• Rao, A. V. "A Survey of Numerical Methods for Optimal Control". Advances in the Astronautical Sciences, volume 135, 2009, pages 497 to 528. AAS paper 09-334. https://www.anilvrao.com/Publications/ConferencePublications/trajectorySurveyAAS.pdf

Vehicle data:

• SpaceX. "Falcon User's Guide". Revision of September 2021. Thrusts, fairing diameter, total vehicle mass and payload capability. https://snebulos.mit.edu/projects/reference/commercial/falcon-users-guide-2021-09.pdf
• Kyle, E. "SpaceX Falcon 9 v1.2 Data Sheet". Space Launch Report, 2018. Stage gross and empty masses and engine specific impulses. These are third-party estimates, not published data, and are treated as such throughout. https://web.archive.org/web/2021/http://www.spacelaunchreport.com/falcon9ft.html

Dependencies:

• numpy 2.0 or later. Trace storage and array reductions. BSD 3-Clause. Harris, C. R. et al. "Array programming with NumPy". Nature, volume 585, 2020, pages 357 to 362. DOI 10.1038/s41586-020-2649-2
• scipy 1.14 or later. scipy.optimize.minimize for SLSQP and Nelder-Mead. BSD 3-Clause. Virtanen, P. et al. "SciPy 1.0: fundamental algorithms for scientific computing in Python". Nature Methods, volume 17, 2020, pages 261 to 272. DOI 10.1038/s41592-019-0686-2
• matplotlib 3.9 or later. Figures, Agg backend only. Matplotlib license, a BSD-compatible license derived from the Python Software Foundation license. Hunter, J. D. "Matplotlib: a 2D graphics environment". Computing in Science and Engineering, volume 9, number 3, 2007, pages 90 to 95. DOI 10.1109/MCSE.2007.55
• pytest 8.3 or later and pytest-cov 6.0 or later. Test runner and coverage measurement. MIT.
• ruff 0.8 or later. Linting and import order. MIT.
• mypy 1.13 or later. Static type checking in strict mode. The package ships a py.typed marker, so an installed copy delivers those types to whatever imports it. MIT.

License

Released under the MIT license. See LICENSE.