Shuttle Reentry#

Mathematical Formulation#

Problem Statement#

Find the optimal angle of attack \(\alpha(t)\) and bank angle \(\beta(t)\) that maximize the crossrange:

(1)#\[J = -\theta(t_f)\]

Subject to the atmospheric entry dynamics:

\[\frac{dh}{dt} = v \sin\gamma\]
\[\frac{d\phi}{dt} = \frac{v \cos\gamma \sin\psi}{r \cos\theta}\]
\[\frac{d\theta}{dt} = \frac{v \cos\gamma \cos\psi}{r}\]
\[\frac{dv}{dt} = -\frac{D}{m} - \frac{\mu \sin\gamma}{r^2}\]
\[\frac{d\gamma}{dt} = \frac{L \cos\beta}{mv} + \cos\gamma\left(\frac{v}{r} - \frac{\mu}{vr^2}\right)\]
\[\frac{d\psi}{dt} = \frac{L \sin\beta}{mv \cos\gamma} + \frac{v \cos\gamma \sin\psi \sin\theta}{r \cos\theta}\]

Where the aerodynamic forces are:

\[L = \frac{1}{2} C_L \rho v^2 S\]
\[D = \frac{1}{2} C_D \rho v^2 S\]
\[C_L = A_0 + A_1 \alpha_{\deg}\]
\[C_D = B_0 + B_1 \alpha_{\deg} + B_2 \alpha_{\deg}^2\]

And the atmospheric model:

\[\rho = \rho_0 e^{-h/H_r}\]
\[\mu = \frac{\mu_E}{r^2}\]

Boundary Conditions#

  • Initial conditions: \(h(0) = 260\) km, \(\phi(0) = 0°\), \(\theta(0) = 0°\), \(v(0) = 25.6\) km/s, \(\gamma(0) = -1°\), \(\psi(0) = 90°\)

  • Final conditions: \(h(t_f) = 80\) km, \(v(t_f) = 2.5\) km/s, \(\gamma(t_f) = -5°\)

  • Control bounds: \(-90° \leq \alpha \leq 90°\), \(-90° \leq \beta \leq 1°\)

  • State bounds: \(h \geq 0\), \(-89° \leq \theta, \gamma, \psi \leq 89°\)

Physical Parameters#

  • Earth gravitational parameter: \(\mu_E = 0.14076539 \times 10^{17}\) ft³/s²

  • Earth radius: \(R_E = 20902900\) ft

  • Reference area: \(S = 2690\) ft²

  • Atmospheric density at sea level: \(\rho_0 = 0.002378\) slug/ft³

  • Density scale height: \(H_r = 23800\) ft

  • Vehicle mass: \(m = 203000/32.174\) slug

  • Aerodynamic coefficients: \(A_0 = -0.20704\), \(A_1 = 0.029244\), \(B_0 = 0.07854\), \(B_1 = -0.61592 \times 10^{-2}\), \(B_2 = 0.621408 \times 10^{-3}\)

State Variables#

  • \(h(t)\): Altitude (ft)

  • \(\phi(t)\): Longitude (rad)

  • \(\theta(t)\): Latitude (rad)

  • \(v(t)\): Velocity magnitude (ft/s)

  • \(\gamma(t)\): Flight path angle (rad)

  • \(\psi(t)\): Azimuth angle (rad)

Control Variables#

  • \(\alpha(t)\): Angle of attack (rad)

  • \(\beta(t)\): Bank angle (rad)

Notes#

This problem optimizes the Space Shuttle’s atmospheric reentry trajectory to maximize crossrange capability while satisfying thermal and structural constraints.

Running This Example#

cd examples/shuttle_reentry
python shuttle_reentry.py

Code Implementation#

examples/shuttle_reentry/shuttle_reentry.py#
  1import casadi as ca
  2import numpy as np
  3
  4import maptor as mtor
  5
  6
  7# ============================================================================
  8# Constants and Scaling
  9# ============================================================================
 10
 11# Physical constants
 12DEG2RAD = np.pi / 180.0
 13MU_EARTH = 0.14076539e17
 14R_EARTH = 20902900.0
 15S_REF = 2690.0
 16RHO0 = 0.002378
 17H_R = 23800.0
 18MASS = 203000.0 / 32.174
 19A0, A1 = -0.20704, 0.029244
 20B0, B1, B2 = 0.07854, -0.61592e-2, 0.621408e-3
 21
 22# Scaling factors for numerical conditioning
 23H_SCALE = 1e5  # Altitude scaling
 24V_SCALE = 1e4  # Velocity scaling
 25
 26# ============================================================================
 27# Problem Setup
 28# ============================================================================
 29
 30problem = mtor.Problem("Shuttle Reentry")
 31phase = problem.set_phase(1)
 32
 33# ============================================================================
 34# Variables
 35# ============================================================================
 36
 37t = phase.time(initial=0.0)
 38h_s = phase.state("altitude_scaled", initial=2.6, final=0.8, boundary=(0, None))
 39phi = phase.state("longitude", initial=0.0)
 40theta = phase.state("latitude", initial=0.0, boundary=(-89 * DEG2RAD, 89 * DEG2RAD))
 41v_s = phase.state("velocity_scaled", initial=2.56, final=0.25, boundary=(1e-4, None))
 42gamma = phase.state(
 43    "flight_path_angle",
 44    initial=-1 * DEG2RAD,
 45    final=-5 * DEG2RAD,
 46    boundary=(-89 * DEG2RAD, 89 * DEG2RAD),
 47)
 48psi = phase.state("heading_angle", initial=90 * DEG2RAD)
 49alpha = phase.control("angle_of_attack", boundary=(-90 * DEG2RAD, 90 * DEG2RAD))
 50beta = phase.control("bank_angle", boundary=(-90 * DEG2RAD, 1 * DEG2RAD))
 51
 52# ============================================================================
 53# Dynamics
 54# ============================================================================
 55
 56# Physics calculations inline for exact numerical behavior
 57h = h_s * H_SCALE
 58v = v_s * V_SCALE
 59r = R_EARTH + h
 60rho = RHO0 * ca.exp(-h / H_R)
 61g = MU_EARTH / r**2
 62alpha_deg = alpha * 180.0 / np.pi
 63CL = A0 + A1 * alpha_deg
 64CD = B0 + B1 * alpha_deg + B2 * alpha_deg**2
 65q = 0.5 * rho * v**2
 66L = q * CL * S_REF
 67D = q * CD * S_REF
 68eps = 1e-10
 69
 70# Dynamics with manual scaling
 71phase.dynamics(
 72    {
 73        h_s: (v * ca.sin(gamma)) / H_SCALE,
 74        phi: (v / r) * ca.cos(gamma) * ca.sin(psi) / (ca.cos(theta) + eps),
 75        theta: (v / r) * ca.cos(gamma) * ca.cos(psi),
 76        v_s: (-(D / MASS) - g * ca.sin(gamma)) / V_SCALE,
 77        gamma: (L / (MASS * v + eps)) * ca.cos(beta) + ca.cos(gamma) * (v / r - g / (v + eps)),
 78        psi: (L * ca.sin(beta) / (MASS * v * ca.cos(gamma) + eps))
 79        + (v / (r * (ca.cos(theta) + eps))) * ca.cos(gamma) * ca.sin(psi) * ca.sin(theta),
 80    }
 81)
 82
 83# ============================================================================
 84# Objective
 85# ============================================================================
 86
 87# Maximize crossrange (final latitude)
 88problem.minimize(-theta.final)
 89
 90# ============================================================================
 91# Mesh and Guess
 92# ============================================================================
 93
 94phase.mesh([8] * 3, np.linspace(-1.0, 1.0, 4))
 95
 96# Initial guess
 97states_guess = []
 98controls_guess = []
 99for N in [8] * 3:
100    t_norm = np.linspace(0, 1, N + 1)
101    h_traj = 2.6 + (0.8 - 2.6) * t_norm
102    phi_traj = np.zeros(N + 1)
103    theta_traj = np.zeros(N + 1)
104    v_traj = 2.56 + (0.25 - 2.56) * t_norm
105    gamma_traj = -1 * DEG2RAD + (-5 * DEG2RAD - (-1 * DEG2RAD)) * t_norm
106    psi_traj = 90 * DEG2RAD * np.ones(N + 1)
107
108    # State order MUST match variable declaration: h_s, phi, theta, v_s, gamma, psi
109    states_guess.append(np.vstack([h_traj, phi_traj, theta_traj, v_traj, gamma_traj, psi_traj]))
110    controls_guess.append(np.vstack([np.zeros(N), -45 * DEG2RAD * np.ones(N)]))
111
112phase.guess(
113    states=states_guess,
114    controls=controls_guess,
115    terminal_time=2000.0,
116)
117
118# ============================================================================
119# Solve
120# ============================================================================
121
122solution = mtor.solve_adaptive(
123    problem,
124    error_tolerance=1e-6,
125    max_iterations=30,
126    min_polynomial_degree=3,
127    max_polynomial_degree=8,
128    nlp_options={
129        "ipopt.max_iter": 2000,
130        "ipopt.tol": 1e-6,
131        "ipopt.constr_viol_tol": 1e-6,
132        "ipopt.acceptable_tol": 1e-3,
133        "ipopt.linear_solver": "mumps",
134        "ipopt.print_level": 0,
135    },
136)
137
138# ============================================================================
139# Results
140# ============================================================================
141
142if solution.status["success"]:
143    crossrange_deg = -solution.status["objective"] * 180.0 / np.pi
144    print(f"Final time: {solution.phases[1]['times']['final']:.1f} seconds")
145    print(f"Crossrange: {crossrange_deg:.2f} degrees")
146    solution.plot()
147else:
148    print(f"Failed: {solution.status['message']}")