Multidisciplinary design optimisation of a small two-stage launcher, taught on a toy problem. Trajectory and propulsion are optimised together with direct collocation (OpenMDAO + Dymos), propellant pairs are compared on equal footing, and a Pareto front is drawn between payload and engine thrust.
Provenance. This repository is inspired by design-optimisation work I did as an R&D engineer in a launch-vehicle start-up. It is a re-implementation of the method on a fictional vehicle: every number here is either a textbook value (Sutton & Biblarz, Rocket Propulsion Elements; NASA CEA; the 1976 Standard Atmosphere) or a generic order of magnitude. Nothing describes a real engine, vehicle, mission or company data.
| Notebook | What you learn | Runtime |
|---|---|---|
00_staging_and_sizing |
Rocket equation and Δv budget, where Isp comes from (CEA tables, ideal-rocket check), a closed-form sizing loop that computes the structural index (hoop stress, shell buckling), a two-objective Pareto front (take-off mass vs length) with NSGA-II | 1 min |
01_two_stage_trajectory_mdao |
The core method: equations of motion per unit frontal area, differentiable engine model from CEA tables, two-phase direct collocation with Dymos, stage separation as a linkage constraint, analytic orbit-insertion closure, validation by explicit re-simulation, mass budget for a chosen diameter | 1 min |
02_propellant_comparison_pareto |
LOX/CH4, LOX/LH2, LOX/RP-1, LOX/HDPE (hybrid) and optionally LOX/HDPE+Al, each flown on its own optimal trajectory; ε-constraint Pareto front payload ↔ maximum thrust with warm starts | 2 min |
The Python package launcher_design/ holds the models; the notebooks only call it.
launcher_design/
atmosphere.py US Standard Atmosphere 1976, self-contained
propellants.py public catalogue of propellant pairs (densities, CEA cards)
thermochem.py CEA look-up tables → scipy interpolants and OpenMDAO metamodels
ideal.py ideal rocket relations (c*, C_F, Isp)
staging.py rocket equation, structural indices, optimal Δv split
sizing.py closed-form two-stage sizing loop + NSGA-II front (pymoo)
ode.py Dymos equations of motion (planar, polar, per unit area)
trajectory.py two-stage MDAO problem: phases, constraints, separation, closure
postprocess.py time series, mass budget for a diameter, dashboard
data/cea_tables/ Ivac, c*, Tc, pc/pe on a (O/F, p_c, ε) grid for each pair (CSV + NPZ)
scripts/ generate_cea_tables.py (needs rocketcea), run_mdao.py (CLI)
tests/ unit tests of the closed-form parts
Write the vehicle per square metre of frontal area: the mass becomes an areal density, thrust and drag become pressures, and the diameter drops out of the trajectory problem (it is chosen in post-processing). Replace the expensive chemical-equilibrium code by tables interpolated with analytic derivatives. Integrate only the boost up to an injection altitude with direct collocation (Gauss-Lobatto, Dymos), and close the orbit insertion analytically (Hohmann + rocket equation): the problem stays small and well conditioned. Optimise the pitch program, the mixture ratios, the nozzle expansion ratios, the engine sizes and the stage split simultaneously under path constraints (acceleration, dynamic pressure) and boundary constraints (lift-off thrust-to-weight, nozzle flow separation at sea level, injection conditions). Then sweep what you cannot put in the objective (thrust cap, propellant pair) to get fronts and comparisons.
Pure Python, no compiler needed (the CEA tables are pre-generated):
python -m venv .venv && source .venv/bin/activate
pip install -e ".[notebooks,dev]"
pytest # closed-form parts
jupyter lab notebooks/ # open 00, 01, 02 in order
python scripts/run_mdao.py LOX_CH4 250000 # CLI: pair, thrust cap [Pa per m²]The optimiser is SciPy's SLSQP through OpenMDAO, so a problem of this size (two phases, 8 segments, ~250 variables) converges in a few seconds on a laptop (tolerance 2e-4 on scaled quantities; tighter tolerances make SLSQP stall on the collocation defects without changing the design, which the explicit re-simulation confirms). Serious work would use SNOPT or IPOPT through pyOptSparse; the code accepts any OpenMDAO driver.
To regenerate the tables (only if you change the propellant catalogue), rocketcea and a Fortran
compiler are needed:
python -m venv .venv-cea && .venv-cea/bin/pip install numpy rocketcea
.venv-cea/bin/python scripts/generate_cea_tables.py- Hybrid pairs (LOX/HDPE, LOX/HDPE+Al) are treated at the thermochemical level only, with the mixture ratio as a free parameter. Regression rate, grain geometry and O/F shift are not modelled.
- Structural indices are inputs of the trajectory problem (0.10 and 0.12 by default) and outputs of the closed-form sizing notebook; coupling the two is left as an exercise.
- No engine-mass model: chamber pressure is therefore fixed, and the upper-stage nozzle goes to its bound. Notebook 02 explains why and what a real study adds.
-
Atmosphere is the 1976 standard, re-implemented (no dependency); drag uses a constant
$C_x$ . - Earth rotation is ignored (non-rotating planet), which costs ~0.4 km/s compared with an eastward launch from a low latitude.
- G. P. Sutton, O. Biblarz, Rocket Propulsion Elements, 9th ed., Wiley, 2017.
- S. Gordon, B. J. McBride, Computer Program for Calculation of Complex Chemical Equilibrium
Compositions and Applications (CEA), NASA RP-1311, 1994. Python wrapper:
rocketcea. - U.S. Standard Atmosphere, 1976, NASA-TM-X-74335.
- R. Falck et al., dymos: A Python package for optimal control of multidisciplinary systems, JOSS 6(59), 2021. J. S. Gray et al., OpenMDAO, Struct. Multidisc. Optim. 59, 2019.
- J. Blank, K. Deb, pymoo: Multi-objective optimization in Python, IEEE Access 8, 2020.
- NASA SP-8007, Buckling of thin-walled circular cylinders, 1968.
MIT.