Space shuttle

Description of the problem

The Space Shuttle reentry problem is a classical benchmark in aerospace optimal control. It models the atmospheric descent of the space shuttle from high altitude to the Terminal Area Energy Management (TAEM) interface. The system includes six state variables: altitude $h(t)$, longitude $\phi(t)$, latitude $\theta(t)$, velocity $v(t)$, flight path angle $\gamma(t)$, and azimuth $\psi(t)$. The control variables are the angle of attack $\alpha(t)$ and bank angle $\beta(t)$. The goal is to maximise the final latitude (crossrange) at TAEM while satisfying aerodynamic, gravitational, and operational constraints [Betts 2010; Bulirsch 1971; Dickmanns 1972].

Mathematical formulation

The problem can be stated as

\[\begin{aligned} \min_{h, \phi, \theta, v, \gamma, \psi, \alpha, \beta} \quad & J(h, \theta) = - \theta(T) \\[1em] \text{s.t.} \quad & \dot{h}(t) = v \sin(\gamma), \\[0.5em] & \dot{\phi}(t) = \frac{v}{r} \cos(\gamma) \frac{\sin(\psi)}{\cos(\theta)}, \\[0.5em] & \dot{\theta}(t) = \frac{v}{r} \cos(\gamma) \cos(\psi), \\[0.5em] & \dot{v}(t) = -\frac{D(h,v,\alpha)}{m} - g(h) \sin(\gamma), \\[0.5em] & \dot{\gamma}(t) = \frac{L(h,v,\alpha)}{m v} \cos(\beta) + \cos(\gamma) \left( \frac{v}{r} - \frac{g(h)}{v} \right), \\[0.5em] & \dot{\psi}(t) = \frac{L(h,v,\alpha)}{m v \cos(\gamma)} \sin(\beta) + \frac{v}{r \cos(\theta)} \cos(\gamma) \sin(\psi) \sin(\theta), \\[0.5em] & \alpha_{\min} \le \alpha(t) \le \alpha_{\max}, \\[0.5em] & \beta_{\min} \le \beta(t) \le \beta_{\max}, \\[0.5em] & h_{\min} \le h(t) \le h_{\max}, \quad v_{\min} \le v(t) \le v_{\max}, \\[0.5em] & \gamma_{\min} \le \gamma(t) \le \gamma_{\max}, \quad \theta_{\min} \le \theta(t) \le \theta_{\max}, \\[0.5em] & h(0) = h_s, \; \phi(0) = \phi_s, \; \theta(0) = \theta_s, \\[0.5em] & v(0) = v_s, \; \gamma(0) = \gamma_s, \; \psi(0) = \psi_s, \\[0.5em] & h(T) = h_t, \; v(T) = v_t, \; \gamma(T) = \gamma_t. \end{aligned}\]

The final time $T$ is free but bounded.

Parameters

Qualitative behaviour

• The optimal trajectory balances lift and drag to control heating and deceleration while extending crossrange.
• The bank angle $\beta$ determines heading changes and crossrange capability.
• The solution typically combines steep reentry to dissipate energy and crossrange manoeuvres to reach the target latitude.

Characteristics

• Nonlinear, six–state dynamics with two controls.
• Strongly nonlinear aerodynamic coefficients.
• Free final time with bounded range.
• Path constraints on states and controls.
• Widely used as a benchmark in optimal control and trajectory optimisation.

References and relevance

• Betts, J.T. (2010). Practical Methods for Optimal Control and Estimation Using Nonlinear Programming. SIAM. Provides a comprehensive discussion of trajectory optimisation for aerospace vehicles, including space shuttle reentry problems. It details numerical methods applicable to nonlinear, constrained, free-final-time problems.

• Bulirsch, R. (1971). Numerical solution of optimal control problems with state constraints by direct methods. Numerische Mathematik. Introduces early direct methods for optimal control problems with constraints, which are foundational for solving shuttle reentry trajectory problems with bounded states and controls.

• Dickmanns, E.D. (1972). Numerical solution methods for nonlinear optimal control problems with state constraints. Automatica. Focuses on numerical techniques for nonlinear optimal control with state constraints, directly relevant to handling the shuttle’s aerodynamic and flight path restrictions.

• Ascher, U.M., Mattheij, R.M.M., & Russell, R.D. (1988). Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. Provides practical algorithms for solving boundary value problems, which underpin the solution of two-point boundary value problems like the shuttle reentry trajectory with fixed initial and terminal states.

• Betts, J.T. (2001). Survey of numerical methods for trajectory optimization. Journal of Guidance, Control, and Dynamics, 24(4), 643–653. Reviews trajectory optimisation methods including direct transcription and collocation, which are commonly applied to the space shuttle reentry benchmark.

Numerical set-up

In this section, we prepare the numerical environment required to study the problem. We begin by importing the relevant Julia packages and then initialise the data frames that will store the results of our simulations and computations. These structures provide the foundation for solving the problem and for comparing the different solution strategies in a consistent way.

using OptimalControlProblems    # to access the Beam model
using OptimalControl            # to import the OptimalControl model
using NLPModelsIpopt            # to solve the model with Ipopt
import DataFrames: DataFrame    # to store data
using NLPModels                 # to retrieve data from the NLP solution
using Plots                     # to plot the trajectories
using Plots.PlotMeasures        # for leftmargin, bottommargin
using JuMP                      # to import the JuMP model
using Ipopt                     # to solve the JuMP model with Ipopt
using Printf                    # to print

data_pb = DataFrame(            # to store data about the problem
    Problem=Symbol[],
    Grid_Size=Int[],
    Variables=Int[],
    Constraints=Int[],
)

data_re = DataFrame(            # to store data about the resolutions
    Model=Symbol[],
    Flag=Any[],
    Iterations=Int[],
    Objective=Float64[],
)

Metadata

The default number of time steps is:

metadata(:space_shuttle)[:grid_size]

500

The default values of the parameters are:

metadata(:space_shuttle)[:parameters]

Parameter = Value
------------------
    t0 =  0.0000e+00
     w =  2.0300e+05
    g₀ =  3.2174e+01
    ρ₀ =  2.3780e-03
    hᵣ =  2.3800e+04
    Rₑ =  2.0903e+07
     μ =  1.4077e+16
     S =  2.6900e+03
    a₀ = -2.0704e-01
    a₁ =  2.9244e-02
    b₀ =  7.8540e-02
    b₁ = -6.1592e-03
    b₂ =  6.2141e-04
Δt_min =  3.5000e+00
Δt_max =  4.5000e+00
  h_t0 =  2.6000e+00
  ϕ_t0 =  0.0000e+00
  θ_t0 =  0.0000e+00
  v_t0 =  2.5600e+00
  γ_t0 = -1.7453e-02
  ψ_t0 =  1.5708e+00
  α_t0 =  0.0000e+00
   α_s =  0.0000e+00
   β_s =  0.0000e+00
  h_tf =  8.0000e-01
  v_tf =  2.5000e-01
  γ_tf = -8.7266e-02
   h_l =  0.0000e+00
   ϕ_l = -6.2832e+00
   ϕ_u =  6.2832e+00
   θ_l = -1.5533e+00
   θ_u =  1.5533e+00
   v_l =  1.0000e-04
   γ_l = -1.5533e+00
   γ_u =  1.5533e+00
   ψ_l = -6.2832e+00
   ψ_u =  6.2832e+00
   α_l = -1.5708e+00
   α_u =  1.5708e+00
   β_l = -1.5533e+00
   β_u =  1.7453e-02

Initial guess

Before solving the problem, it is often useful to inspect the initial guess (sometimes called the first iterate). This guess is obtained by running the NLP solver with max_iter = 0, which evaluates the problem formulation without performing any optimisation steps.

We plot the resulting trajectories for both the OptimalControl and JuMP models. Since both backends represent the same mathematical problem, their initial guesses should coincide, providing a useful consistency check before moving on to the optimised solution.

function plot_initial_guess(problem)

    # -----------------------------
    # Build OptimalControl problem
    # -----------------------------
    docp = eval(problem)(OptimalControlBackend())
    nlp_oc = nlp_model(docp)
    ocp_oc = ocp_model(docp)

    # Solve NLP with zero iterations (initial guess)
    nlp_oc_sol = NLPModelsIpopt.ipopt(nlp_oc; max_iter=0)

    # Build OptimalControl solution
    ocp_sol = build_ocp_solution(docp, nlp_oc_sol)

    # get dimensions
    n = state_dimension(ocp_oc)
    m = control_dimension(ocp_oc)

    # -----------------------------
    # Plot OptimalControl solution
    # -----------------------------
    plt = plot(
        ocp_sol;
        state_style=(color=1,),
        costate_style=(color=1, legend=:none),
        control_style=(color=1, legend=:none),
        path_style=(color=1, legend=:none),
        dual_style=(color=1, legend=:none),
        size=(816, 220*(n+m)),
        label="OptimalControl",
        leftmargin=20mm,
    )

    # Hide legend for additional state plots
    for i in 2:n
        plot!(plt[i]; legend=:none)
    end

    # -----------------------------
    # Build JuMP model
    # -----------------------------
    nlp_jp = eval(problem)(JuMPBackend())

    # Solve NLP with zero iterations (initial guess)
    set_optimizer(nlp_jp, Ipopt.Optimizer)
    set_optimizer_attribute(nlp_jp, "max_iter", 0)
    optimize!(nlp_jp)

    # Extract trajectories
    t_grid = time_grid(nlp_jp)
    x_fun = state(nlp_jp)
    u_fun = control(nlp_jp)
    p_fun = costate(nlp_jp)

    # -----------------------------
    # Plot JuMP solution on top
    # -----------------------------
    # States
    for i in 1:n
        label = i == 1 ? "JuMP" : :none
        plot!(plt[i], t_grid, t -> x_fun(t)[i]; color=2, linestyle=:dash, label=label)
    end

    # Costates
    for i in 1:n
        plot!(plt[n+i], t_grid, t -> -p_fun(t)[i]; color=2, linestyle=:dash, label=:none)
    end

    # Controls
    for i in 1:m
        plot!(plt[2*n+i], t_grid, t -> u_fun(t)[i]; color=2, linestyle=:dash, label=:none)
    end

    return plt
end

plot_initial_guess(:space_shuttle)

Solving the problem

To solve an optimal control problem, we can rely on two complementary formulations: the OptimalControl backend, which works directly with the discretised control problem, and the JuMP backend, which leverages JuMP’s flexible modelling framework.

Both approaches generate equivalent NLPs that can be solved with Ipopt, and comparing them ensures consistency between the two formulations.

Before solving, we can inspect the discretisation details of the problem. The table below reports the number of grid points, decision variables, and constraints associated with the chosen formulation.

push!(data_pb,(
    Problem=:space_shuttle,
    Grid_Size=metadata(:space_shuttle)[:grid_size],
    Variables=get_nvar(nlp_model(space_shuttle(OptimalControlBackend()))),
    Constraints=get_ncon(nlp_model(space_shuttle(OptimalControlBackend()))),
))

OptimalControl model

We first solve the problem using the OptimalControl backend. The process begins by importing the problem definition and constructing the associated nonlinear programming (NLP) model. This NLP is then passed to the Ipopt solver, with standard options for tolerance and barrier parameter strategy.

# import DOCP model
docp = space_shuttle(OptimalControlBackend())

# get NLP model
nlp_oc = nlp_model(docp)

# solve
nlp_oc_sol = NLPModelsIpopt.ipopt(
    nlp_oc;
    print_level=4,
    tol=1e-8,
    mu_strategy="adaptive",
    sb="yes",
)

Total number of variables............................:     4009
                     variables with only lower bounds:     1002
                variables with lower and upper bounds:     3007
                     variables with only upper bounds:        0
Total number of equality constraints.................:     3009
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


Number of Iterations....: 109

                                   (scaled)                 (unscaled)
Objective...............:  -5.9587501338421633e-01    5.9587501338421633e-01
Dual infeasibility......:   1.3122876897850429e-10    1.3122876897850429e-10
Constraint violation....:   2.6527414709320851e-09    2.6527414709320851e-09
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   5.7619865070546199e-11    5.7619865070546199e-11
Overall NLP error.......:   2.6527414709320851e-09    2.6527414709320851e-09


Number of objective function evaluations             = 110
Number of objective gradient evaluations             = 110
Number of equality constraint evaluations            = 110
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 110
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 109
Total seconds in IPOPT                               = 5.695

EXIT: Optimal Solution Found.

JuMP model

We now repeat the procedure using the JuMP backend. Here, the problem is reformulated as a JuMP model, which offers a flexible and widely used framework for nonlinear optimisation in Julia. The solver settings are chosen to mirror those used previously, so that the results can be compared on an equal footing.

# import model
nlp_jp = space_shuttle(JuMPBackend())

# solve with Ipopt
set_optimizer(nlp_jp, Ipopt.Optimizer)
set_optimizer_attribute(nlp_jp, "print_level", 4)
set_optimizer_attribute(nlp_jp, "tol", 1e-8)
set_optimizer_attribute(nlp_jp, "mu_strategy", "adaptive")
set_optimizer_attribute(nlp_jp, "linear_solver", "mumps")
set_optimizer_attribute(nlp_jp, "sb", "yes")
optimize!(nlp_jp)

Total number of variables............................:     4009
                     variables with only lower bounds:     1002
                variables with lower and upper bounds:     3007
                     variables with only upper bounds:        0
Total number of equality constraints.................:     3009
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


Number of Iterations....: 115

                                   (scaled)                 (unscaled)
Objective...............:  -5.9587501342513571e-01    5.9587501342513571e-01
Dual infeasibility......:   3.4751664545807872e-10    3.4751664545807872e-10
Constraint violation....:   7.0468874763252032e-09    7.0468874763252032e-09
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.3610483608806734e-10    1.3610483608806734e-10
Overall NLP error.......:   7.0468874763252032e-09    7.0468874763252032e-09


Number of objective function evaluations             = 116
Number of objective gradient evaluations             = 116
Number of equality constraint evaluations            = 116
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 116
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 115
Total seconds in IPOPT                               = 7.328

EXIT: Optimal Solution Found.

Numerical comparisons

In this section, we examine the results of the problem resolutions. We extract the solver status (flag), the number of iterations, and the objective value for each model. This provides a first overview of how each approach performs and sets the stage for a more detailed comparison of the solution trajectories.

# from OptimalControl model
push!(data_re,(
    Model=:OptimalControl,
    Flag=nlp_oc_sol.status,
    Iterations=nlp_oc_sol.iter,
    Objective=nlp_oc_sol.objective,
))

# from JuMP model
push!(data_re,(
    Model=:JuMP,
    Flag=termination_status(nlp_jp),
    Iterations=barrier_iterations(nlp_jp),
    Objective=objective_value(nlp_jp),
))

We compare the solutions obtained from the OptimalControl and JuMP models by examining the number of iterations required for convergence, the $L^2$-norms of the differences in states, controls, and additional variables, and the corresponding objective values. Both absolute and relative errors are reported, providing a clear quantitative measure of the agreement between the two approaches.

function L2_norm(T, X)
    # T and X are supposed to be one dimensional
    s = 0.0
    for i in 1:(length(T) - 1)
        s += 0.5 * (X[i]^2 + X[i + 1]^2) * (T[i + 1]-T[i])
    end
    return √(s)
end

function print_numerical_comparisons(problem, docp, nlp_oc_sol, nlp_jp)

    # get relevant data from OptimalControl model
    ocp_sol = build_ocp_solution(docp, nlp_oc_sol)
    t_oc = time_grid(ocp_sol)
    x_oc = state(ocp_sol).(t_oc)
    u_oc = control(ocp_sol).(t_oc)
    v_oc = variable(ocp_sol)
    o_oc = objective(ocp_sol)
    i_oc = iterations(ocp_sol)

    # get relevant data from JuMP model
    t_jp = time_grid(nlp_jp)
    x_jp = state(nlp_jp).(t_jp)
    u_jp = control(nlp_jp).(t_jp)
    o_jp = objective(nlp_jp)
    v_jp = variable(nlp_jp)
    i_jp = iterations(nlp_jp)

    x_vars = state_components(nlp_jp)
    u_vars = control_components(nlp_jp)
    v_vars = variable_components(nlp_jp)

    println("┌─ ", string(problem))
    println("│")
    println("├─  Number of Iterations")
    @printf("│     OptimalControl : %d   JuMP : %d\n", i_oc, i_jp)

    # States
    println("├─  States (L2 Norms)")
    for i in eachindex(x_vars)
        xi_oc = [x_oc[k][i] for k in eachindex(t_oc)]
        xi_jp = [x_jp[k][i] for k in eachindex(t_jp)]
        L2_ae = L2_norm(t_oc, xi_oc - xi_jp)
        L2_re = L2_ae / (0.5 * (L2_norm(t_oc, xi_oc) + L2_norm(t_oc, xi_jp)))
        @printf("│     %-6s Abs: %.3e   Rel: %.3e\n", x_vars[i], L2_ae, L2_re)
    end

    # Controls
    println("├─  Controls (L2 Norms)")
    for i in eachindex(u_vars)
        ui_oc = [u_oc[k][i] for k in eachindex(t_oc)]
        ui_jp = [u_jp[k][i] for k in eachindex(t_jp)]
        L2_ae = L2_norm(t_oc, ui_oc - ui_jp)
        L2_re = L2_ae / (0.5 * (L2_norm(t_oc, ui_oc) + L2_norm(t_oc, ui_jp)))
        @printf("│     %-6s Abs: %.3e   Rel: %.3e\n", u_vars[i], L2_ae, L2_re)
    end

    # Variables
    if !isnothing(v_vars)
        println("├─  Variables")
        for i in eachindex(v_vars)
            vi_oc = v_oc[i]
            vi_jp = v_jp[i]
            vi_ae = abs(vi_oc - vi_jp)
            vi_re = vi_ae / (0.5 * (abs(vi_oc) + abs(vi_jp)))
            @printf("│     %-6s Abs: %.3e   Rel: %.3e\n", v_vars[i], vi_ae, vi_re)
        end
    end

    # Objective
    o_ae = abs(o_oc - o_jp)
    o_re = o_ae / (0.5 * (abs(o_oc) + abs(o_jp)))
    println("├─  Objective")
    @printf("│            Abs: %.3e   Rel: %.3e\n", o_ae, o_re)
    println("└─")
    return nothing
end

print_numerical_comparisons(:space_shuttle, docp, nlp_oc_sol, nlp_jp)

┌─ space_shuttle
│
├─  Number of Iterations
│     OptimalControl : 109   JuMP : 115
├─  States (L2 Norms)
│     scaled_h Abs: 1.894e-08   Rel: 2.393e-10
│     ϕ      Abs: 8.215e-09   Rel: 1.889e-10
│     θ      Abs: 2.783e-09   Rel: 1.860e-10
│     scaled_v Abs: 1.268e-08   Rel: 1.631e-10
│     γ      Abs: 2.494e-08   Rel: 2.730e-08
│     ψ      Abs: 5.974e-07   Rel: 1.287e-08
├─  Controls (L2 Norms)
│     α      Abs: 4.630e-08   Rel: 3.403e-09
│     β      Abs: 2.170e-05   Rel: 6.036e-07
├─  Variables
│     tf     Abs: 4.877e-07   Rel: 2.428e-10
├─  Objective
│            Abs: 4.092e-11   Rel: 6.867e-11
└─

Plotting the solutions

In this section, we visualise the trajectories of the states, costates, and controls obtained from both the OptimalControl and JuMP solutions. The plots provide an intuitive way to compare the two approaches and to observe how the constraints and the optimal control influence the system dynamics.

For each variable, the OptimalControl solution is shown in solid lines, while the JuMP solution is overlaid using dashed lines. Since both models represent the same mathematical problem, their trajectories should closely coincide, highlighting the consistency between the two formulations.

# build an ocp solution to use the plot from OptimalControl package
ocp_sol = build_ocp_solution(docp, nlp_oc_sol)

# dimensions
n = state_dimension(ocp_sol)
m = control_dimension(ocp_sol)

# from OptimalControl solution
plt = plot(
    ocp_sol;
    color=1,
    size=(816, 240*(n+m)),
    label="OptimalControl",
    leftmargin=20mm,
)
for i in 2:length(plt)
    plot!(plt[i]; legend=:none)
end

# from JuMP solution
t = time_grid(nlp_jp)     # t0, ..., tN = tf
x = state(nlp_jp)         # function of time
u = control(nlp_jp)       # function of time
p = costate(nlp_jp)       # function of time

for i in 1:n # state
    label = i == 1 ? "JuMP" : :none
    plot!(plt[i], t, t -> x(t)[i]; color=2, linestyle=:dash, label=label)
end

for i in 1:n # costate
    plot!(plt[n+i], t, t -> -p(t)[i]; color=2, linestyle=:dash, label=:none)
end

for i in 1:m # control
    plot!(plt[2n+i], t, t -> u(t)[i]; color=2, linestyle=:dash, label=:none)
end