Maximizing Payload Fraction Under the Tsiolkovsky Equation
Introduction
Every rocket designer eventually runs into the same brutal arithmetic: the Tsiolkovsky rocket equation punishes you exponentially for wanting more velocity change. Add propellant, and you also add structure to hold that propellant, which in turn demands even more propellant to move. Multi-stage rockets exist precisely to escape this trap — by discarding dead structural mass along the way, you avoid dragging spent tanks and engines all the way to orbit.
But staging introduces a new question: if a mission requires a fixed total velocity change (Δv), how should that Δv be split among the stages? Give too much to a low-efficiency booster and you waste propellant; give too much to a high-efficiency upper stage and you carry excess structure through the atmosphere unnecessarily. This is the optimal staging problem, and today we’ll solve it concretely: for a fixed target Δv, we want the Δv allocation across stages that maximizes the payload fraction (equivalently, minimizes propellant and structural mass for a given payload).
Mathematical Formulation
The Tsiolkovsky equation
For a single stage, the velocity change is:
$$
\Delta v = I_{sp} , g_0 , \ln\left(\frac{m_0}{m_f}\right)
$$
where $I_{sp}$ is the specific impulse, $g_0$ is standard gravity, and $m_0/m_f$ is the mass ratio $R$.
Stage payload fraction
Each stage $i$ has a structural coefficient $\sigma_i$, the fraction of that stage’s own hardware+propellant mass that is dead weight (tanks, engines, avionics) rather than propellant:
$$
\sigma_i = \frac{m_{struct,i}}{m_{struct,i} + m_{prop,i}}
$$
If $\pi_i = m_{PL,i}/m_{0,i}$ is the fraction of stage $i$’s total initial mass that is “payload” (everything riding above it, including later stages), then the mass ratio and payload fraction are linked by:
$$
R_i = \frac{1}{\sigma_i(1-\pi_i) + \pi_i}
\quad\Longrightarrow\quad
\pi_i = \frac{\dfrac{1}{R_i} - \sigma_i}{1 - \sigma_i}
$$
Because each stage’s payload is the assembly of every stage above it, the fractions telescope: the fraction of the total liftoff mass that reaches orbit as pure payload is simply the product across all stages:
$$
\Pi = \prod_{i=1}^{N} \pi_i(\Delta v_i)
$$
The optimization problem
$$
\max_{\Delta v_1, \dots, \Delta v_N} ; \Pi = \prod_{i=1}^{N} \pi_i(\Delta v_i)
\qquad \text{subject to} \qquad
\sum_{i=1}^{N} \Delta v_i = \Delta v_{total}
$$
Taking logs turns the product into a sum, which is far friendlier for numerical solvers:
$$
\max ; \sum_{i=1}^{N} \ln \pi_i(\Delta v_i)
$$
At the optimum, Lagrange’s condition tells us the marginal effectiveness of adding one more m/s of Δv must be equal across every stage:
$$
\frac{\partial \ln \pi_i}{\partial \Delta v_i} = \lambda \quad \text{(same } \lambda \text{ for all } i\text{)}
$$
This is the numerical fingerprint we’ll check to confirm the solver actually found the true optimum.
A Concrete Example: 3-Stage Orbital Launcher
We target a total Δv of 9500 m/s — a realistic figure for reaching low Earth orbit once gravity and drag losses are folded in. Our three stages use progressively more efficient (and more expensive) propulsion:
| Stage | Propellant type | $I_{sp}$ [s] | $\sigma$ (structural coefficient) |
|---|---|---|---|
| 1 | Kerolox booster | 282 | 0.06 |
| 2 | Kerolox second stage | 311 | 0.08 |
| 3 | Hydrolox upper stage | 450 | 0.11 |
Intuitively, the high-$I_{sp}$ hydrolox stage should be “worth more” Δv than the lower-$I_{sp}$ booster, but its higher structural coefficient pulls in the opposite direction. The optimizer resolves this trade-off numerically.
Python Source Code (Google Colaboratory)
1 | # ============================================================ |
Code Walkthrough
Section 1 — Constants and specs. We fix standard gravity, the total mission Δv, and pack each stage’s $I_{sp}$ and $\sigma$ into small NumPy arrays. Keeping these as arrays (rather than separate variables) lets every downstream function operate on all three stages at once without loops.
Section 2 — Physics core. mass_ratio is a direct implementation of the Tsiolkovsky equation solved for $R = e^{\Delta v/(I_{sp} g_0)}$. stage_payload_fraction inverts the relationship derived above to recover $\pi_i$ from a chosen $\Delta v_i$. total_payload_fraction multiplies the per-stage fractions together, implementing the telescoping product $\Pi = \prod \pi_i$.
Section 3 — The optimizer. Rather than maximizing the product directly (numerically unstable when payload fractions are small), we minimize the negative sum of logs, which is mathematically equivalent but far more numerically stable. SLSQP (Sequential Least Squares Programming) is the right choice here because it natively supports both bounds and an equality constraint — exactly what “$\sum \Delta v_i = \Delta v_{total}$” requires. The penalty branch (return 1e6) protects the optimizer from ever evaluating the objective at an infeasible point where a stage’s payload fraction would go negative, which would otherwise throw a math domain error from log.
Section 4 — Optimality verification. This is the empirical check of the Lagrange condition. We nudge each $\Delta v_i$ by ±1 m/s and take a central difference of $\ln \pi_i$. If the optimizer truly found the constrained optimum, all three marginal effectiveness values $\lambda_i$ should converge to (nearly) the same number — proof that no further improvement is possible by shifting Δv between any pair of stages.
Section 5 — Reporting. Plain print statements summarize the optimal allocation, the resulting payload fraction, and the percentage improvement over a naive equal three-way split of Δv.
Section 6 — Visualization. The 3D surface sweeps $\Delta v_1$ and $\Delta v_2$ over a grid (with $\Delta v_3$ determined by the constraint), computing the resulting payload fraction at every grid point — fully vectorized with NumPy broadcasting rather than nested Python loops, so the entire 70×70 grid evaluates essentially instantly. The optimum found by SLSQP is plotted as a highlighted marker sitting on the peak of the surface, visually confirming the numerical result. The bar chart puts the optimal and equal-split Δv allocations side by side for a direct, human-readable comparison.
Runtime Performance
This problem has only three decision variables and a smooth, well-behaved objective, so SLSQP converges in well under a second — no acceleration techniques are needed for the optimization itself. The only part that touches any real volume of computation is the 3D surface grid, and vectorizing it with NumPy array operations (as done above, with zero Python-level loops over grid points) keeps that well under a second as well.
============================================================ OPTIMAL STAGING RESULT ============================================================ Stage 1: dv = 1110.04 m/s | pi_1 = 0.648284 | marginal d(ln pi)/d(dv) = -0.00039720 Stage 2: dv = 2378.88 m/s | pi_2 = 0.411314 | marginal d(ln pi)/d(dv) = -0.00039720 Stage 3: dv = 6011.08 m/s | pi_3 = 0.164173 | marginal d(ln pi)/d(dv) = -0.00039720 ------------------------------------------------------------ Total delta-v check : 9500.000 m/s (target 9500.0) Optimal payload fraction : 4.3776 % Equal-split payload frac.: 3.4746 % Improvement over equal split: 25.99 % ============================================================
Interpreting the Graphs
The 3D surface (left panel) shows the payload fraction as a smooth hill over the $(\Delta v_1, \Delta v_2)$ plane, with $\Delta v_3$ implicitly filling in whatever Δv remains. The surface falls off sharply near the edges — where one stage is forced to carry almost all of the mission’s Δv and its mass ratio blows up exponentially — and rises to a single interior peak. The cyan marker sits exactly on that peak, confirming that SLSQP located the true maximum rather than a boundary artifact or local irregularity.
The bar chart (right panel) translates the abstract optimum into an engineering decision: how many m/s of Δv should each real stage actually be designed to deliver? Because the hydrolox third stage has the highest $I_{sp}$, expect the optimizer to shift more of the mission’s total Δv onto it relative to a naive equal split — the exponential penalty for low-$I_{sp}$ mass ratios makes the lower stages relatively more “expensive” per unit of Δv, even though their structural coefficients are smaller.
Conclusion
The optimal staging problem is a clean example of constrained nonlinear optimization hiding inside a deceptively simple physical formula. By reformulating the payload-fraction product as a log-sum, applying SLSQP under an equality constraint, and cross-checking the result against the analytic Lagrange condition, we get both a numerically verified optimum and an intuitive visualization of the entire solution landscape. The same framework generalizes directly to four or more stages, asymmetric Δv-loss profiles (gravity and drag losses concentrated in the lower stages), or reusable first-stage penalties — all by extending the ISP and SIGMA arrays and letting the optimizer do the rest.
