Aero-Assisted Plane Change#

Mathematical Formulation#

Problem Statement#

Find the optimal lift coefficient \(C_L(t)\) and bank angle \(\beta(t)\) that maximize the final velocity:

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

Subject to the atmospheric flight dynamics:

\[\frac{d\phi}{dt} = \frac{v \cos\gamma \sin\psi}{r}\]
\[\frac{dh}{dt} = v \sin\gamma\]
\[\frac{dv}{dt} = -a_1 \rho v^2 (1 + C_L^2) - \frac{\mu \sin\gamma}{r^2}\]
\[\frac{d\gamma}{dt} = a_0 \rho v (C_L \cos\beta + M \cos\gamma)\]
\[\frac{d\psi}{dt} = \frac{a_0 \rho v C_L \sin\beta}{\cos\gamma} - \frac{v \cos\gamma \cos\psi \tan\phi}{r}\]

Where:

\[r = R_e + h\]
\[\rho = \rho_0 e^{-(h-h_0)/h_r}\]
\[v_s = \sqrt{\frac{\mu}{R_e}}\]
\[\rho_s = \rho_0 e^{h_0/h_r}\]
\[M = \frac{1}{a_0 \rho r}\left(1 - \frac{\mu}{rv^2}\right)\]
\[\dot{q} = 17600 \sqrt{\frac{\rho}{\rho_s}} \left(\frac{v}{v_s}\right)^{3.15}\]

Heat Rate Constraint#

\[\dot{q}(t) \leq 800 \text{ BTU/ft}^2\text{-s}\]

Terminal Constraint#

\[\cos\phi(t_f) \cos\psi(t_f) = \cos(18°)\]

Boundary Conditions#

  • Initial conditions: \(\phi(0) = 0°\), \(h(0) = 365000\) ft, \(v(0) = 25745.704\) ft/s, \(\gamma(0) = -0.55°\), \(\psi(0) = 0°\)

  • Final conditions: \(h(t_f) = 365000\) ft, others determined by optimization

  • Control bounds: \(0 \leq C_L \leq 2\), \(0 \leq \beta \leq 180°\)

  • State bounds: \(h \geq 0\), \(20000 \leq v \leq 28000\) ft/s, \(-89° \leq \phi, \gamma, \psi \leq 89°\)

Physical Parameters#

  • Earth radius: \(R_e = 2.092643 \times 10^7\) ft

  • Vehicle mass: \(m = 3.315 \times 10^2\) slug

  • Sea level density: \(\rho_0 = 3.3195 \times 10^{-5}\) slug/ft³

  • Reference altitude: \(h_0 = 10^5\) ft

  • Scale height: \(h_r = 2.41388 \times 10^4\) ft

  • Drag coefficient: \(C_{D0} = 0.032\)

  • Lift-drag parameter: \(k = 1.4\)

  • Reference area: \(S = 1.2584 \times 10^2\) ft²

  • Gravitational parameter: \(\mu = 1.40895 \times 10^{16}\) ft³/s²

  • Derived parameters: \(a_0 = \frac{S}{2m}\sqrt{\frac{C_{D0}}{k}}\), \(a_1 = \frac{C_{D0} S}{2m}\)

State Variables#

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

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

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

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

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

Control Variables#

  • \(C_L(t)\): Lift coefficient

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

Notes#

This problem involves using atmospheric forces to change orbital inclination while managing heat rate constraints and maximizing final velocity.

Running This Example#

cd examples/aero-assisted_plane_change
python aero-assisted_plane_change.py

Code Implementation#

examples/aero-assisted_plane_change/aero-assisted_plane_change.py#
  1import casadi as ca
  2import numpy as np
  3
  4import maptor as mtor
  5
  6
  7# Constants from Table 10.3
  8Re = 2.092643e7  # ft
  9m = 3.315e2  # slugs
 10rho0 = 3.3195e-5  # slug/ft³
 11h0 = 1e5  # ft
 12hr = 2.41388e4  # ft
 13CD0 = 0.032
 14k = 1.4
 15S = 1.2584e2  # ft²
 16mu = 1.40895e16  # ft³/sec²
 17
 18DEG2RAD = np.pi / 180.0
 19
 20# Derived parameters
 21rhos = rho0 * ca.exp(h0 / hr)
 22a0 = (S / (2 * m)) * ca.sqrt(CD0 / k)
 23a1 = (CD0 * S) / (2 * m)
 24vs = ca.sqrt(mu / Re)
 25
 26# Problem setup
 27problem = mtor.Problem("Optimal Aero-Assisted Plane Change")
 28phase = problem.set_phase(1)
 29
 30# Variables
 31t = phase.time(initial=0.0, final=(800, 2000))
 32phi = phase.state(
 33    "phi", initial=0.0, final=(-89 * DEG2RAD, 89 * DEG2RAD), boundary=(-89 * DEG2RAD, 89 * DEG2RAD)
 34)
 35h = phase.state("altitude", initial=365000.0, final=365000.0, boundary=(0, 400000))
 36v = phase.state("velocity", initial=25745.704, boundary=(20000, 28000))
 37gamma = phase.state(
 38    "gamma",
 39    initial=-0.55 * DEG2RAD,
 40    final=(-10 * DEG2RAD, 10 * DEG2RAD),
 41    boundary=(-10 * DEG2RAD, 10 * DEG2RAD),
 42)
 43psi = phase.state(
 44    "psi", initial=0.0, final=(-89 * DEG2RAD, 89 * DEG2RAD), boundary=(-89 * DEG2RAD, 89 * DEG2RAD)
 45)
 46
 47CL = phase.control("lift_coeff", boundary=(0, 2))
 48beta = phase.control("bank_angle", boundary=(0, 180 * DEG2RAD))
 49
 50# Additional variables
 51r = Re + h
 52rho = rho0 * ca.exp(-(h - h0) / hr)
 53M = (1 / (a0 * rho * r)) * (1 - mu / (r * v**2))
 54q_dot = 17600 * ca.sqrt(rho / rhos) * (v / vs) ** 3.15
 55
 56# Dynamics
 57phi_dot = (v / r) * ca.cos(gamma) * ca.sin(psi)
 58h_dot = v * ca.sin(gamma)
 59v_dot = -a1 * rho * v**2 * (1 + CL**2) - (mu * ca.sin(gamma)) / r**2
 60gamma_dot = a0 * rho * v * (CL * ca.cos(beta) + M * ca.cos(gamma))
 61psi_dot = (a0 * rho * v * CL * ca.sin(beta)) / ca.cos(gamma) - (
 62    v * ca.cos(gamma) * ca.cos(psi) * ca.tan(phi)
 63) / r
 64
 65phase.dynamics(
 66    {
 67        phi: phi_dot,
 68        h: h_dot,
 69        v: v_dot,
 70        gamma: gamma_dot,
 71        psi: psi_dot,
 72    }
 73)
 74
 75# Path constraints: heat rate limits (applied at every collocation point)
 76phase.path_constraints(
 77    q_dot <= 800,  # Heat rate upper limit
 78    q_dot >= 0,  # Heat rate lower limit (physical constraint)
 79)
 80
 81# Event constraint: final boundary condition (applied at final boundary only)
 82target_angle = ca.cos(18 * DEG2RAD)
 83phase.event_constraints(ca.cos(phi.final) * ca.cos(psi.final) == target_angle)
 84
 85# Objective: Maximize final velocity
 86problem.minimize(-v.final)
 87
 88# Mesh and guess
 89phase.mesh([12, 12, 12], [-1.0, -1 / 3, 1 / 3, 1.0])
 90
 91states_guess = []
 92controls_guess = []
 93for N in [12, 12, 12]:
 94    tau = np.linspace(-1, 1, N + 1)
 95    t_norm = (tau + 1) / 2
 96
 97    # Linear interpolation for states
 98    phi_vals = 0.0 + (0.0 - 0.0) * t_norm  # phi stays near 0
 99    h_vals = 365000.0 + (365000.0 - 365000.0) * t_norm  # altitude constant
100    v_vals = 25715.704 + (25715.704 - 25715.704) * t_norm  # velocity constant
101    gamma_vals = -0.55 * DEG2RAD + (0.0 - (-0.55 * DEG2RAD)) * t_norm
102    psi_vals = 0.0 + (18 * DEG2RAD - 0.0) * t_norm  # approach target angle
103
104    states_guess.append(np.vstack([phi_vals, h_vals, v_vals, gamma_vals, psi_vals]))
105
106    # Control guess
107    CL_vals = np.full(N, 1.0)  # Mid-range lift coefficient
108    beta_vals = np.full(N, 45 * DEG2RAD)  # Mid-range bank angle
109    controls_guess.append(np.vstack([CL_vals, beta_vals]))
110
111phase.guess(
112    states=states_guess,
113    controls=controls_guess,
114    terminal_time=800,
115)
116
117# Solve
118solution = mtor.solve_adaptive(
119    problem,
120    error_tolerance=1e-6,
121    max_iterations=50,
122    min_polynomial_degree=4,
123    max_polynomial_degree=15,
124    nlp_options={
125        "ipopt.print_level": 5,
126        "ipopt.max_iter": 3000,
127        "ipopt.tol": 1e-8,
128        "ipopt.constr_viol_tol": 1e-7,
129    },
130)
131
132# Results
133if solution.status["success"]:
134    final_velocity = -solution.status["objective"]  # Convert back from minimization
135    print(f"Final velocity: {final_velocity:.4f} ft/sec")
136    print(f"Final time: {solution.phases[1]['times']['final']:.4f} sec")
137    print("Reference: v_f = 22043.5079 ft/sec, t_f = 1005.8778 sec")
138
139    error_v = abs(final_velocity - 22043.5079) / 22043.5079 * 100
140    print(f"Velocity error: {error_v:.2f}%")
141
142    solution.plot()
143else:
144    print(f"Failed: {solution.status['message']}")