---
title: Space Shuttle Reentry Trajectory
---

**Originally Contributed by**: Henrique Ferrolho

This tutorial demonstrates how to compute a reentry trajectory for the
[Space Shuttle](https://en.wikipedia.org/wiki/Space_Shuttle),
by formulating and solving a nonlinear programming problem.
The problem was drawn from Chapter 6 of 
["Practical Methods for Optimal Control and Estimation Using Nonlinear Programming"](https://epubs.siam.org/doi/book/10.1137/1.9780898718577),
by John T. Betts.

The motion of the vehicle is defined by the following set of DAEs:

$$
\begin{align}
\dot{h} & = v \sin \gamma , \\
\dot{\phi} & = \frac{v}{r} \cos \gamma \sin \psi / \cos \theta , \\
\dot{\theta} & = \frac{v}{r} \cos \gamma \cos \psi , \\
\dot{v} & = -\frac{D}{m} - g \sin \gamma , \\
\dot{\gamma} & = \frac{L}{m v} \cos(\beta) + \cos \gamma \left ( \frac{v}{r} - \frac{g}{v} \right ) , \\
\dot{\psi} & = \frac{1}{m v \cos \gamma} L \sin(\beta) + \frac{v}{r \cos \theta} \cos \gamma \sin \psi \sin \theta , \\
q & \le q_U , \\
\end{align}
$$

where the aerodynamic heating on the vehicle wing leading edge is $q = q_a q_r$
and the dynamic variables are

$$
\begin{align}
     h & \quad \text{altitude (ft)},     \qquad & & \gamma \quad \text{flight path angle (rad)}, \\
  \phi & \quad \text{longitude (rad)},   \qquad & & \psi   \quad \text{azimuth (rad)},           \\
\theta & \quad \text{latitude (rad)},    \qquad & & \alpha \quad \text{angle of attack (rad)},   \\
     v & \quad \text{velocity (ft/sec)}, \qquad & & \beta  \quad \text{bank angle (rad)}.
\end{align}
$$

The aerodynamic and atmospheric forces on the vehicle are
specified by the following quantities (English units):

$$
\begin{align}
           D & = \frac{1}{2} c_D S \rho v^2,                  \qquad & a_0 & = -0.20704, \\
           L & = \frac{1}{2} c_L S \rho v^2,                  \qquad & a_1 & =  0.029244, \\
           g & = \mu / r^2,                                   \qquad & \mu & =  0.14076539 \times 10^{17}, \\
           r & = R_e + h,                                     \qquad & b_0 & =  0.07854, \\
        \rho & = \rho_0 \exp[-h/h_r],                         \qquad & b_1 & = -0.61592  \times 10^{-2}, \\
      \rho_0 & = 0.002378,                                    \qquad & b_2 & =  0.621408 \times 10^{-3}, \\
         h_r & = 23800,                                       \qquad & q_r & =  17700 \sqrt{\rho} (0.0001 v)^{3.07}, \\
         c_L & = a_0 + a_1 \hat{\alpha},                      \qquad & q_a & =  c_0 + c_1 \hat{\alpha} + c_2 \hat{\alpha}^2 + c_3 \hat{\alpha}^3, \\
         c_D & = b_0 + b_1 \hat{\alpha} + b_2 \hat{\alpha}^2, \qquad & c_0 & =  1.0672181, \\
\hat{\alpha} & = 180 \alpha / \pi,                            \qquad & c_1 & = -0.19213774 \times 10^{-1}, \\
         R_e & = 20902900,                                    \qquad & c_2 & =  0.21286289 \times 10^{-3}, \\
           S & = 2690,                                        \qquad & c_3 & = -0.10117249 \times 10^{-5}.
\end{align}
$$

The reentry trajectory begins at an altitude where the aerodynamic forces are quite small
with a weight of $w = 203000$ (lb) and mass $m = w / g_0$ (slug), where $g_0 = 32.174$ (ft/sec$^2$).
The initial conditions are as follows:

$$
\begin{align}
     h & = 260000 \text{ ft},  \qquad & v      & = 25600 \text{ ft/sec}, \\
  \phi & = 0      \text{ deg}, \qquad & \gamma & = -1    \text{ deg},    \\
\theta & = 0      \text{ deg}, \qquad & \psi   & = 90    \text{ deg}.
\end{align}
$$

The final point on the reentry trajectory occurs at the unknown (free) time $t_F$,
at the so-called terminal area energy management (TAEM) interface, which is defined by the conditions

$$
h = 80000 \text{ ft}, \qquad v = 2500 \text{ ft/sec}, \qquad \gamma = -5 \text{ deg}.
$$

As explained in the book, our goal is to maximize the final crossrange,
which is equivalent to maximizing the final latitude of the vehicle,
i.e., $J = \theta(t_F)$.



### Approach 
We will use a discretized model of time, with a fixed number of discretized points, $n$.
The decision variables at each point are going to be the state of the vehicle and the controls commanded to it.
In addition, we will also make each time step size $\Delta t$ a decision variable;
that way, we can either fix the time step size easily,
or allow the solver to fine-tune the duration between each adjacent pair of points.
Finally, in order to approximate the derivatives of the problem dynamics, we will use either rectangular
or [trapezoidal](http://en.wikipedia.org/wiki/Trapezoidal_rule) integration.

### Disclaimer
Do not try to actually land a Space Shuttle using this notebook! 😛
There's no mesh refinement going on, which can lead to unrealistic trajectories having position
and velocity errors with orders of magnitude $10^4$ ft and $10^2$ ft/sec, respectively.


```julia
using Interpolations, Plots
using JuMP, Ipopt  # , KNITRO
```


```julia
# Global variables
const w  = 203000.0  # weight (lb)
const g₀ = 32.174    # acceleration (ft/sec^2)
const m  = w / g₀    # mass (slug)

# Aerodynamic and atmospheric forces on the vehicle
const ρ₀ =  0.002378
const hᵣ =  23800.0
const Rₑ =  20902900.0
const μ  =  0.14076539e17
const S  =  2690.0
const a₀ = -0.20704
const a₁ =  0.029244
const b₀ =  0.07854
const b₁ = -0.61592e-2
const b₂ =  0.621408e-3
const c₀ =  1.0672181
const c₁ = -0.19213774e-1
const c₂ =  0.21286289e-3
const c₃ = -0.10117249e-5

# Initial conditions
const hₛ = 2.6          # altitude (ft) / 1e5
const ϕₛ = deg2rad(0)   # longitude (rad)
const θₛ = deg2rad(0)   # latitude (rad)
const vₛ = 2.56         # velocity (ft/sec) / 1e4
const γₛ = deg2rad(-1)  # flight path angle (rad)
const ψₛ = deg2rad(90)  # azimuth (rad)
const αₛ = deg2rad(0)   # angle of attack (rad)
const βₛ = deg2rad(0)   # banck angle (rad)
const tₛ = 1.00         # time step (sec)

# Final conditions, the so-called Terminal Area Energy Management (TAEM)
const hₜ = 0.8          # altitude (ft) / 1e5
const vₜ = 0.25         # velocity (ft/sec) / 1e4
const γₜ = deg2rad(-5)  # flight path angle (rad)

# Number of mesh points (knots) to be used
const n = 2009

# Integration scheme to be used for the dynamics
const integration_rule = "rectangular"
```




    "rectangular"




```julia
user_options_ipopt = (
    "mu_strategy" => "monotone",
    "linear_solver" => "ma27",  # For the best results, it is advised to experiment different linear solvers.
                                # If Ipopt is not compiled with MA27/MA57, it may fallback to 'MUMPS'.
                                # In general, the linear solver MA27 is much faster than MUMPS.
)

# user_options_knitro = (
#     "algorithm" => KNITRO.KN_ALG_BAR_DIRECT,
#     "bar_murule" => KNITRO.KN_BAR_MURULE_QUALITY,
#     "linsolver" => KNITRO.KN_LINSOLVER_MA27,
# )

# Create JuMP model, using Ipopt as the solver
model = Model(optimizer_with_attributes(Ipopt.Optimizer, user_options_ipopt...))
# model = Model(optimizer_with_attributes(KNITRO.Optimizer, user_options_knitro...))

@variables(model, begin
               0 ≤ scaled_h[1:n]                # altitude (ft) / 1e5
                          ϕ[1:n]                # longitude (rad)
    deg2rad(-89) ≤        θ[1:n] ≤ deg2rad(89)  # latitude (rad)
            1e-4 ≤ scaled_v[1:n]                # velocity (ft/sec) / 1e4
    deg2rad(-89) ≤        γ[1:n] ≤ deg2rad(89)  # flight path angle (rad)
                          ψ[1:n]                # azimuth (rad)
    deg2rad(-90) ≤        α[1:n] ≤ deg2rad(90)  # angle of attack (rad)
    deg2rad(-89) ≤        β[1:n] ≤ deg2rad( 1)  # banck angle (rad)
    #        0.5 ≤       Δt[1:n] ≤ 1.5          # time step (sec)
                         Δt[1:n] == 1.0         # time step (sec)
end)

# Fix initial conditions
fix(scaled_h[1], hₛ; force=true)
fix(       ϕ[1], ϕₛ; force=true)
fix(       θ[1], θₛ; force=true)
fix(scaled_v[1], vₛ; force=true)
fix(       γ[1], γₛ; force=true)
fix(       ψ[1], ψₛ; force=true)

# Fix final conditions
fix(scaled_h[n], hₜ; force=true)
fix(scaled_v[n], vₜ; force=true)
fix(       γ[n], γₜ; force=true)

# Initial guess: linear interpolation between boundary conditions
xₛ = [hₛ, ϕₛ, θₛ, vₛ, γₛ, ψₛ, αₛ, βₛ, tₛ]
xₜ = [hₜ, ϕₛ, θₛ, vₜ, γₜ, ψₛ, αₛ, βₛ, tₛ]
interp_linear = LinearInterpolation([1, n], [xₛ, xₜ])
initial_guess = mapreduce(transpose, vcat, interp_linear.(1:n))
set_start_value.(all_variables(model), vec(initial_guess))

# Functions to restore `h` and `v` to their true scale
@NLexpression(model, h[j = 1:n], scaled_h[j] * 1e5)
@NLexpression(model, v[j = 1:n], scaled_v[j] * 1e4)

# Helper functions
# @NLexpression(model, c_L[j = 1:n], a₀ + a₁ * rad2deg(α[j]))
# @NLexpression(model, c_D[j = 1:n], b₀ + b₁ * rad2deg(α[j]) + b₂ * rad2deg(α[j])^2)
@NLexpression(model, c_L[j = 1:n], a₀ + a₁ * (180 * α[j] / π))
@NLexpression(model, c_D[j = 1:n], b₀ + b₁ * (180 * α[j] / π) + b₂ * (180 * α[j] / π)^2)
@NLexpression(model,   ρ[j = 1:n], ρ₀ * exp(-h[j] / hᵣ))
@NLexpression(model,   D[j = 1:n], 0.5 * c_D[j] * S * ρ[j] * v[j]^2)
@NLexpression(model,   L[j = 1:n], 0.5 * c_L[j] * S * ρ[j] * v[j]^2)
@NLexpression(model,   r[j = 1:n], Rₑ + h[j])
@NLexpression(model,   g[j = 1:n], μ / r[j]^2)

# Motion of the vehicle as a differential-algebraic system of equations (DAEs)
@NLexpression(model, δh[j = 1:n], v[j] * sin(γ[j]))
@NLexpression(model, δϕ[j = 1:n], (v[j] / r[j]) * cos(γ[j]) * sin(ψ[j]) / cos(θ[j]))
@NLexpression(model, δθ[j = 1:n], (v[j] / r[j]) * cos(γ[j]) * cos(ψ[j]))
@NLexpression(model, δv[j = 1:n], -(D[j] / m) - g[j] * sin(γ[j]))
@NLexpression(model, δγ[j = 1:n], (L[j] / (m * v[j])) * cos(β[j]) + cos(γ[j]) * ((v[j] / r[j]) - (g[j] / v[j])))
@NLexpression(model, δψ[j = 1:n], (1 / (m * v[j] * cos(γ[j]))) * L[j] * sin(β[j]) + (v[j] / (r[j] * cos(θ[j]))) * cos(γ[j]) * sin(ψ[j]) * sin(θ[j]))

# System dynamics
for j in 2:n
    i = j - 1  # index of previous knot

    if integration_rule == "rectangular"
        # Rectangular integration
        @NLconstraint(model, h[j] == h[i] + Δt[i] * δh[i])
        @NLconstraint(model, ϕ[j] == ϕ[i] + Δt[i] * δϕ[i])
        @NLconstraint(model, θ[j] == θ[i] + Δt[i] * δθ[i])
        @NLconstraint(model, v[j] == v[i] + Δt[i] * δv[i])
        @NLconstraint(model, γ[j] == γ[i] + Δt[i] * δγ[i])
        @NLconstraint(model, ψ[j] == ψ[i] + Δt[i] * δψ[i])
    elseif integration_rule == "trapezoidal"
        # Trapezoidal integration
        @NLconstraint(model, h[j] == h[i] + 0.5 * Δt[i] * (δh[j] + δh[i]))
        @NLconstraint(model, ϕ[j] == ϕ[i] + 0.5 * Δt[i] * (δϕ[j] + δϕ[i]))
        @NLconstraint(model, θ[j] == θ[i] + 0.5 * Δt[i] * (δθ[j] + δθ[i]))
        @NLconstraint(model, v[j] == v[i] + 0.5 * Δt[i] * (δv[j] + δv[i]))
        @NLconstraint(model, γ[j] == γ[i] + 0.5 * Δt[i] * (δγ[j] + δγ[i]))
        @NLconstraint(model, ψ[j] == ψ[i] + 0.5 * Δt[i] * (δψ[j] + δψ[i]))
    else
        @error "Unexpected integration rule '$(integration_rule)'"
    end
end

# Objective: Maximize crossrange
@objective(model, Max, θ[n])

# Solve for the control and state
status = optimize!(model)

# Show final crossrange of the solution
println("Final latitude θ = ", round(objective_value(model) |> rad2deg, digits = 2), "°")
```

    
    ******************************************************************************
    This program contains Ipopt, a library for large-scale nonlinear optimization.
     Ipopt is released as open source code under the Eclipse Public License (EPL).
             For more information visit http://projects.coin-or.org/Ipopt
    ******************************************************************************
    
    This is Ipopt version 3.13.2, running with linear solver ma27.
    
    Number of nonzeros in equality constraint Jacobian...:    72260
    Number of nonzeros in inequality constraint Jacobian.:        0
    Number of nonzeros in Lagrangian Hessian.............:   162574
    
    Total number of variables............................:    16063
                         variables with only lower bounds:     4014
                    variables with lower and upper bounds:     8033
                         variables with only upper bounds:        0
    Total number of equality constraints.................:    12048
    Total number of inequality constraints...............:        0
            inequality constraints with only lower bounds:        0
       inequality constraints with lower and upper bounds:        0
            inequality constraints with only upper bounds:        0
    
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
       0  0.0000000e+00 6.54e+02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
       1 -6.3343098e-02 5.67e+03 1.69e+03  -1.0 1.83e+00    -  6.10e-01 1.00e+00f  1
       2 -6.9086603e-02 2.07e+03 3.21e+04  -1.0 7.82e-01   4.0 4.04e-01 1.00e+00h  1
       3 -4.7521502e-02 1.06e+03 3.25e+04  -1.0 6.32e-01   4.4 7.74e-01 6.88e-01h  1
       4 -4.7196157e-02 1.04e+03 3.23e+04  -1.0 4.67e-01   4.9 1.00e+00 1.40e-02h  1
       5 -4.6646809e-02 1.02e+03 3.26e+04  -1.0 4.99e-01   5.3 1.00e+00 2.28e-02h  1
       6 -4.3199203e-02 9.01e+02 2.88e+04  -1.0 4.65e-01   4.8 1.00e+00 1.21e-01h  1
       7 -4.1589219e-02 8.51e+02 2.86e+04  -1.0 4.56e-01   5.2 1.00e+00 5.69e-02h  1
       8 -3.9815973e-02 7.99e+02 5.04e+04  -1.0 4.61e-01   5.7 1.00e+00 6.26e-02h  1
       9 -3.5837693e-02 6.94e+02 3.29e+04  -1.0 5.45e-01   5.2 1.00e+00 1.37e-01h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      10 -3.4315783e-02 6.54e+02 4.95e+04  -1.0 4.15e-01   5.6 1.00e+00 5.93e-02h  1
      11 -3.2793480e-02 6.14e+02 1.57e+05  -1.0 3.74e-01   6.0 1.00e+00 6.18e-02h  1
      12 -2.9817475e-02 5.43e+02 7.32e+04  -1.0 5.89e-01   5.6 1.00e+00 1.21e-01h  1
      13 -2.8732303e-02 5.17e+02 1.44e+05  -1.0 3.73e-01   6.0 1.00e+00 4.94e-02h  1
      14 -2.7624819e-02 4.91e+02 4.58e+05  -1.0 3.19e-01   6.4 1.00e+00 5.23e-02h  1
      15 -2.5514458e-02 4.44e+02 1.45e+05  -1.0 4.34e-01   5.9 1.00e+00 9.94e-02h  1
      16 -2.4689804e-02 4.25e+02 4.06e+05  -1.0 2.88e-01   6.4 1.00e+00 4.28e-02h  1
      17 -2.2692027e-02 3.84e+02 2.17e+05  -1.0 8.96e-01   5.9 1.00e+00 1.02e-01h  1
      18 -2.1971616e-02 3.68e+02 3.68e+05  -1.0 3.13e-01   6.3 1.00e+00 4.11e-02h  1
      19 -2.1205108e-02 3.52e+02 1.21e+06  -1.0 2.74e-01   6.7 1.00e+00 4.52e-02h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      20 -1.9734474e-02 3.23e+02 3.13e+05  -1.0 3.22e-01   6.2 1.00e+00 8.63e-02h  1
      21 -1.9141687e-02 3.11e+02 1.06e+06  -1.0 2.62e-01   6.7 1.00e+00 3.79e-02h  1
      22 -1.7735591e-02 2.84e+02 3.63e+05  -1.0 3.32e-01   6.2 1.00e+00 8.89e-02h  1
      23 -1.7194352e-02 2.74e+02 9.42e+05  -1.0 2.72e-01   6.6 1.00e+00 3.76e-02h  1
      24 -1.5896638e-02 2.50e+02 4.02e+05  -1.0 3.81e-01   6.1 1.00e+00 8.89e-02h  1
      25 -1.5419422e-02 2.42e+02 8.63e+05  -1.0 2.93e-01   6.6 1.00e+00 3.60e-02h  1
      26 -1.4255760e-02 2.22e+02 4.51e+05  -1.0 6.09e-01   6.1 1.00e+00 8.63e-02h  1
      27 -1.3815887e-02 2.14e+02 8.35e+05  -1.0 3.13e-01   6.5 1.00e+00 3.58e-02h  1
      28 -1.3557670e-02 2.09e+02 6.95e+05  -1.0 2.69e+00   6.0 1.81e-01 2.08e-02h  1
      29 -1.3253507e-02 2.04e+02 8.32e+05  -1.0 3.66e-01   6.5 1.00e+00 2.50e-02h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      30 -1.2564402e-02 1.93e+02 2.93e+06  -1.0 3.33e-01   6.9 1.00e+00 5.79e-02h  1
      31 -1.1900903e-02 1.82e+02 7.18e+05  -1.0 4.22e-01   6.4 1.00e+00 5.54e-02h  1
      32 -1.1555694e-02 1.77e+02 2.99e+06  -1.0 3.56e-01   6.8 1.00e+00 3.04e-02h  1
      33 -1.0947335e-02 1.68e+02 6.60e+05  -1.0 4.58e-01   6.4 1.00e+00 5.27e-02h  1
      34 -1.0544195e-02 1.62e+02 3.13e+06  -1.0 3.66e-01   6.8 1.00e+00 3.67e-02h  1
      35 -1.0095108e-02 1.56e+02 1.03e+06  -1.0 5.35e-01   6.3 1.00e+00 3.89e-02h  1
      36 -9.7865301e-03 1.51e+02 4.68e+06  -1.0 4.95e-01   6.7 1.00e+00 2.76e-02h  1
      37 -9.1324034e-03 1.43e+02 1.27e+06  -1.0 5.66e-01   6.3 1.00e+00 5.70e-02h  1
      38 -8.8501547e-03 1.39e+02 5.15e+06  -1.0 3.95e-01   6.7 1.00e+00 2.60e-02h  1
      39 -8.0375525e-03 1.30e+02 1.13e+06  -1.0 5.82e-01   6.2 1.00e+00 7.34e-02h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      40 -7.6240699e-03 1.25e+02 4.36e+06  -1.0 3.52e-01   6.6 1.00e+00 4.05e-02h  1
      41 -7.1438175e-03 1.19e+02 1.36e+07  -1.0 3.42e-01   7.1 1.00e+00 4.83e-02h  1
      42 -6.1719132e-03 1.08e+02 3.11e+06  -1.0 3.54e-01   6.6 9.41e-01 9.84e-02h  1
      43 -5.3081700e-03 9.77e+01 8.49e+06  -1.0 3.12e-01   7.0 1.00e+00 9.73e-02h  1
      44 -3.6286106e-03 8.02e+01 6.45e+06  -1.0 2.96e-01   6.5 2.38e-01 2.03e-01h  1
      45 -2.9676043e-03 7.26e+01 6.18e+06  -1.0 2.28e-01   7.0 9.41e-04 9.97e-02h  1
      46 -2.7545976e-03 7.01e+01 6.28e+06  -1.0 3.54e-01   6.5 3.61e-04 3.49e-02h  1
      47 -2.7284479e-03 6.98e+01 1.07e+06  -1.0 2.37e-01   6.9 1.00e+00 4.50e-03h  1
      48 -6.3723607e-04 5.07e+01 1.08e+06  -1.0 4.06e-01   6.4 3.35e-02 3.55e-01h  1
      49 -5.3811689e-04 4.94e+01 9.13e+05  -1.0 3.65e-01   6.9 1.00e+00 2.65e-02h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      50  7.5437165e-04 4.58e+01 7.31e+05  -1.0 4.78e-01   6.4 3.67e-01 3.47e-01h  1
      51  9.5641099e-04 4.23e+01 7.94e+05  -1.0 3.66e-01   6.8 1.00e+00 8.85e-02h  1
      52  1.8339791e-03 3.14e+01 6.61e+05  -1.0 3.83e-01   6.3 5.76e-01 4.23e-01h  1
      53  2.5314070e-03 1.18e+01 3.66e+05  -1.0 2.83e-01   5.9 7.30e-02 6.11e-01f  1
      54  2.5566378e-03 1.11e+01 3.46e+05  -1.0 6.55e-01   5.4 9.32e-02 5.92e-02f  1
      55  2.9678499e-03 2.17e+00 1.53e+05  -1.0 1.10e-01   5.8 2.53e-01 1.00e+00f  1
      56  2.9715632e-03 6.26e-03 4.62e+03  -1.0 4.87e-03   5.3 9.85e-01 1.00e+00f  1
      57  3.0331020e-03 4.58e-07 1.39e+02  -1.0 1.96e-03   4.9 1.00e+00 1.00e+00f  1
      58  3.5602871e-03 7.08e-05 1.17e+02  -1.0 5.56e-03   4.4 1.00e+00 8.75e-01f  1
      59  3.7601790e-03 8.87e-06 1.78e+01  -1.0 2.26e-03   3.9 1.00e+00 1.00e+00f  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      60  4.6968224e-03 2.01e-04 1.89e+01  -1.0 7.19e-03   3.4 1.00e+00 1.00e+00f  1
      61  5.6733696e-03 3.15e-04 6.72e+00  -1.0 7.67e-03   2.9 1.00e+00 1.00e+00f  1
      62  8.1993558e-03 3.12e-03 5.67e+00  -1.0 1.94e-02   2.5 1.00e+00 1.00e+00f  1
      63  1.2104241e-02 3.84e-02 2.87e+00  -1.0 2.94e-02   2.0 1.00e+00 1.00e+00f  1
      64  2.0189014e-02 4.28e-01 1.92e+00  -1.0 5.92e-02   1.5 1.00e+00 1.00e+00f  1
      65  2.3315131e-02 5.18e-02 1.97e+00  -1.0 2.27e-02   1.9 1.00e+00 1.00e+00f  1
      66  3.0070243e-02 6.30e-01 1.39e+00  -1.0 4.82e-02   1.5 1.00e+00 1.00e+00f  1
      67  3.2625400e-02 7.66e-02 1.37e+00  -1.0 1.79e-02   1.9 1.00e+00 1.00e+00f  1
      68  3.8859978e-02 9.16e-01 1.27e+00  -1.0 4.93e-02   1.4 1.00e+00 1.00e+00f  1
      69  4.1165937e-02 1.08e-01 1.08e+00  -1.0 1.60e-02   1.8 1.00e+00 1.00e+00f  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      70  4.7087350e-02 1.30e+00 1.28e+00  -1.0 5.59e-02   1.4 1.00e+00 1.00e+00f  1
      71  4.9277379e-02 1.47e-01 1.05e+00  -1.0 1.75e-02   1.8 1.00e+00 1.00e+00f  1
      72  5.5047421e-02 1.78e+00 1.26e+00  -1.0 6.16e-02   1.3 1.00e+00 1.00e+00f  1
      73  5.7189693e-02 1.93e-01 9.86e-01  -1.0 1.85e-02   1.7 1.00e+00 1.00e+00f  1
      74  5.8340406e-02 8.98e-02 2.03e-01  -1.7 1.29e-02   1.3 1.00e+00 1.00e+00f  1
      75  6.1889239e-02 8.08e-01 2.30e-01  -1.7 3.82e-02   0.8 1.00e+00 1.00e+00f  1
      76  7.0954393e-02 8.80e+00 2.31e-01  -1.7 1.20e-01   0.3 1.00e+00 1.00e+00f  1
      77  7.4368103e-02 5.72e-01 1.16e-01  -1.7 2.31e-02   0.7 1.00e+00 1.00e+00f  1
      78  7.5645870e-02 1.58e-01 2.07e-02  -2.5 1.36e-02   0.3 1.00e+00 1.00e+00h  1
      79  7.5958778e-02 5.91e-03 1.34e-03  -3.8 2.56e-03  -0.2 1.00e+00 1.00e+00h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      80  7.6281162e-02 2.80e-03 3.51e-04  -5.7 1.81e-03  -0.7 1.00e+00 1.00e+00h  1
      81  7.7260527e-02 2.59e-02 3.62e-04  -5.7 5.50e-03  -1.2 1.00e+00 1.00e+00h  1
      82  8.0233183e-02 2.54e-01 3.74e-04  -5.7 1.70e-02  -1.7 1.00e+00 1.00e+00h  1
      83  8.9503188e-02 2.93e+00 4.11e-04  -5.7 5.60e-02  -2.1 1.00e+00 1.00e+00h  1
      84  9.3292307e-02 3.30e-01 3.11e-04  -5.7 1.61e-02  -1.7 1.00e+00 1.00e+00h  1
      85  1.0483151e-01 4.67e+00 4.23e-04  -5.7 6.34e-02  -2.2 1.00e+00 1.00e+00h  1
      86  1.0987479e-01 4.95e-01 4.81e-04  -5.7 2.78e-02  -1.8 1.00e+00 1.00e+00h  1
      87  1.3186068e-01 6.43e+02 2.05e-02  -5.7 9.80e-01  -2.2 4.00e-01 1.00e+00h  1
      88  1.4668209e-01 2.42e+02 1.15e-02  -5.7 1.72e-01  -1.8 1.00e+00 1.00e+00h  1
      89  1.5770031e-01 1.77e+02 8.44e-03  -5.7 2.84e-01  -2.3 1.00e+00 2.81e-01h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      90  1.8834000e-01 4.18e+02 7.79e-03  -5.7 1.19e+01  -2.8 2.72e-02 6.33e-02h  1
      91  2.4430759e-01 1.55e+02 1.52e-03  -5.7 2.91e-01  -2.3 1.00e+00 1.00e+00h  1
      92  2.7261698e-01 3.20e+02 1.17e-03  -5.7 1.82e+00  -2.8 3.72e-01 2.59e-01h  1
      93  2.7590990e-01 2.88e+02 5.55e-03  -5.7 2.91e-01  -2.4 1.00e+00 9.83e-02h  1
      94  3.0194885e-01 2.27e+02 4.84e-03  -5.7 5.44e+00  -2.9 5.76e-02 2.21e-01h  1
      95  3.1894809e-01 1.26e+02 3.16e-03  -5.7 2.97e-01  -2.4 4.40e-01 4.52e-01h  1
      96  3.2036792e-01 1.07e+02 1.04e-02  -5.7 2.13e-01  -2.0 1.00e+00 1.54e-01h  1
      97  3.4549672e-01 3.72e+01 9.33e-04  -5.7 3.70e-01  -2.5 9.55e-01 6.83e-01h  1
      98  3.4958894e-01 2.18e+01 2.49e-03  -5.7 1.14e-01  -2.1 1.00e+00 4.11e-01h  1
      99  3.7186722e-01 1.35e+01 8.03e-04  -5.7 3.61e-01  -2.5 5.77e-01 6.79e-01h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
     100  3.8126592e-01 6.84e-01 9.59e-04  -5.7 1.28e-01  -2.1 1.00e+00 1.00e+00h  1
     101  3.9410111e-01 1.73e+00 9.95e-04  -5.7 4.90e-01  -2.6 4.45e-01 4.84e-01h  1
     102  3.9791537e-01 9.97e-01 2.88e-03  -5.7 1.34e-01  -2.2 1.00e+00 4.29e-01h  1
     103  4.2267620e-01 1.20e+01 1.06e-03  -5.7 6.21e-01  -2.6 3.63e-01 1.00e+00h  1
     104  4.2925878e-01 9.36e-01 3.63e-03  -5.7 1.52e-01  -2.2 6.64e-01 1.00e+00h  1
     105  4.3549818e-01 1.93e+00 2.10e-03  -5.7 8.64e-01  -2.7 1.52e-01 2.94e-01h  1
     106  4.4352082e-01 5.80e-01 8.49e-04  -5.7 1.68e-01  -2.3 1.00e+00 1.00e+00f  1
     107  4.5709315e-01 4.57e+01 4.96e-03  -5.7 2.60e+00  -2.7 1.05e-01 4.96e-01h  1
     108  4.6818779e-01 6.20e+00 6.21e-03  -5.7 1.06e+00  -2.3 1.00e+00 1.00e+00h  1
     109  4.9736625e-01 1.35e+01 6.43e-03  -5.7 6.85e-01  -2.8 9.16e-02 1.00e+00h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
     110  5.0752610e-01 2.96e+00 9.90e-04  -5.7 1.19e-01  -2.4 8.24e-01 1.00e+00h  1
     111  5.3306382e-01 9.96e+00 1.05e-03  -5.7 7.88e-01  -2.9 1.84e-01 1.00e+00h  1
     112  5.4009482e-01 2.59e+00 3.86e-04  -5.7 1.16e-01  -2.4 7.55e-01 1.00e+00h  1
     113  5.5492270e-01 1.95e+01 3.25e-04  -5.7 3.61e-01  -2.9 7.67e-01 1.00e+00h  1
     114  5.7867469e-01 1.21e+02 2.16e-03  -5.7 9.40e-01  -3.4 3.69e-01 1.00e+00h  1
     115  6.0105011e-01 2.46e+02 3.34e-02  -5.7 6.42e+00  -3.9 5.65e-02 1.82e-01h  1
     116  6.0876620e-01 6.35e+01 1.54e-02  -5.7 1.00e+00  -2.5 9.71e-01 9.62e-01h  1
     117  6.2419028e-01 4.88e+01 1.07e-02  -5.7 2.62e+00    -  6.29e-02 3.29e-01h  1
     118  6.1074851e-01 6.93e+01 5.39e-03  -5.7 7.24e-01  -3.0 7.27e-02 1.00e+00h  1
     119  5.9969204e-01 1.69e+02 2.96e-03  -5.7 7.49e-01  -3.5 4.65e-01 1.00e+00h  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
     120  5.9578494e-01 8.96e+00 1.19e-03  -5.7 4.27e-01  -4.0 6.04e-01 1.00e+00H  1
     121  5.9650165e-01 9.43e+00 3.72e-04  -5.7 3.00e-01    -  7.87e-01 1.00e+00h  1
     122  5.9599152e-01 3.63e-02 1.20e-05  -5.7 1.18e-01    -  1.00e+00 1.00e+00h  1
     123  5.9598532e-01 1.52e-04 1.75e-06  -5.7 3.13e-02    -  1.00e+00 1.00e+00h  1
     124  5.9607738e-01 9.02e-03 1.49e-06  -8.6 6.56e-02    -  9.59e-01 1.00e+00h  1
     125  5.9604860e-01 5.48e-05 2.72e-07  -8.6 2.80e-02    -  1.00e+00 1.00e+00h  1
     126  5.9604675e-01 2.07e-06 7.25e-08  -8.6 1.48e-02    -  1.00e+00 1.00e+00h  1
     127  5.9604651e-01 2.55e-07 1.27e-08  -8.6 6.20e-03    -  1.00e+00 1.00e+00h  1
     128  5.9604649e-01 1.81e-08 9.01e-10  -8.6 1.65e-03    -  1.00e+00 1.00e+00h  1
     129  5.9604649e-01 2.14e-08 4.59e-10 -12.9 1.18e-03    -  9.99e-01 1.00e+00h  1
    
    Number of Iterations....: 129
    
                                       (scaled)                 (unscaled)
    Objective...............:  -5.9604648822552675e-01    5.9604648822552675e-01
    Dual infeasibility......:   4.5924084974112041e-10    4.5924084974112041e-10
    Constraint violation....:   9.2127296347221943e-09    2.1434971131384373e-08
    Complementarity.........:   1.6911024325149909e-10   -1.6911024325149909e-10
    Overall NLP error.......:   9.2127296347221943e-09    2.1434971131384373e-08
    
    
    Number of objective function evaluations             = 132
    Number of objective gradient evaluations             = 130
    Number of equality constraint evaluations            = 132
    Number of inequality constraint evaluations          = 0
    Number of equality constraint Jacobian evaluations   = 130
    Number of inequality constraint Jacobian evaluations = 0
    Number of Lagrangian Hessian evaluations             = 129
    Total CPU secs in IPOPT (w/o function evaluations)   =     12.277
    Total CPU secs in NLP function evaluations           =     30.827
    
    EXIT: Optimal Solution Found.
    Final latitude θ = 34.15°



```julia
ts = cumsum([0; value.(Δt)])[1:end-1]

plt_altitude = plot(ts, value.(scaled_h), legend = nothing, title = "Altitude (100,000 ft)")
plt_longitude = plot(ts, rad2deg.(value.(ϕ)), legend = nothing, title = "Longitude (deg)")
plt_latitude = plot(ts, rad2deg.(value.(θ)), legend = nothing, title = "Latitude (deg)")
plt_velocity = plot(ts, value.(scaled_v), legend = nothing, title = "Velocity (1000 ft/sec)")
plt_flight_path = plot(ts, rad2deg.(value.(γ)), legend = nothing, title = "Flight Path (deg)")
plt_azimuth = plot(ts, rad2deg.(value.(ψ)), legend = nothing, title = "Azimuth (deg)")

plt = plot(plt_altitude,  plt_velocity, plt_longitude, plt_flight_path, plt_latitude,
           plt_azimuth, layout=grid(3, 2), linewidth=2, size=(700, 700))
```




    

    




```julia
function q(h, v, a)
    ρ(h) = ρ₀ * exp(-h / hᵣ)
    qᵣ(h, v) = 17700 * √ρ(h) * (0.0001 * v)^3.07
    qₐ(a) = c₀ + c₁ * rad2deg(a) + c₂ * rad2deg(a)^2 + c₃ * rad2deg(a)^3    
    # Aerodynamic heating on the vehicle wing leading edge
    qₐ(a) * qᵣ(h, v)
end

plt_attack_angle = plot(ts[1:end-1], rad2deg.(value.(α)[1:end-1]), legend=nothing, title="Angle of Attack (deg)")
plt_bank_angle = plot(ts[1:end-1], rad2deg.(value.(β)[1:end-1]), legend=nothing, title="Bank Angle (deg)")
plt_heating = plot(ts, q.(value.(scaled_h)*1e5, value.(scaled_v)*1e4, value.(α)), legend=nothing, title="Heating (BTU/ft/ft/sec)") 

plt = plot(plt_attack_angle, plt_bank_angle, plt_heating, layout=grid(3, 1), linewidth=2, size=(700, 700))
```




    

    




```julia
plt = plot(rad2deg.(value.(ϕ)), rad2deg.(value.(θ)), value.(scaled_h),
           linewidth=2, legend=nothing, title="Space Shuttle Reentry Trajectory",
           xlabel="Longitude (deg)", ylabel="Latitude (deg)", zlabel="Altitude (100,000 ft)")
```




    

    



4004ce996af70086804ad2fef4f60778a0051b09
ipynb
Jupyter Notebook
notebook/modelling/space_shuttle_reentry_trajectory.ipynb
ferrolho/JuMPTutorials.jl
3e50e6f12d0542fd598d1cdbe0c602a6ae5ea271
MIT
notebook/modelling/space_shuttle_reentry_trajectory.ipynb
ferrolho/JuMPTutorials.jl
3e50e6f12d0542fd598d1cdbe0c602a6ae5ea271
MIT
notebook/modelling/space_shuttle_reentry_trajectory.ipynb
ferrolho/JuMPTutorials.jl
3e50e6f12d0542fd598d1cdbe0c602a6ae5ea271
MIT
Qwen/Qwen-72B
1. YES
2. YES