Goddard with state constraints

This well-known problem^{[1] [2]} models the ascent of a rocket through the atmosphere, and we restrict here ourselves to vertical (one dimensional) trajectories. The state variables are the altitude $r$, speed $v$ and mass $m$ of the rocket during the flight, for a total dimension of 3. The rocket is subject to gravity $g$, thrust $u$ and drag force $D$ (function of speed and altitude). The final time $T$ is free, and the objective is to reach a maximal altitude with a bounded fuel consumption.

We thus want to solve the optimal control problem in Mayer form

\[ \max\, r(T)\]

subject to the control dynamics

\[ \dot{r} = v, \quad \dot{v} = \frac{T_{\max}\,u - D(r,v)}{m} - g, \quad \dot{m} = -u,\]

and subject to the control constraint $u(t) \in [0,1]$ and the state constraint $v(t) \leq v_{\max}$. The initial state is fixed while only the final mass is prescribed.

Model and solution

using CTProblems
using DifferentialEquations
using Plots

You can access the problem in the CTProblems package:

prob = Problem(:goddard)

title           = Goddard problem with state constraint - maximise altitude
model    (Type) = OptimalControlModel{Autonomous, NonFixed}
solution (Type) = OptimalControlSolution

Then, the model is given by

prob.model

The (autonomous) optimal control problem is given by:

    Cd = 310
    Tmax = 3.5
    β = 500
    b = 2
    t0 = 0
    r0 = 1
    v0 = 0
    vmax = 0.1
    m0 = 1
    mf = 0.6
    x0 = [r0, v0, m0]
    tf ∈ R, variable
    t ∈ [t0, tf], time
    x ∈ R³, state
    u ∈ R, control
    r = x₁
    v = x₂
    m = x₃
    0 ≤ u(t) ≤ 1, u_con
    r(t) ≥ r0, x_con_rmin
    0 ≤ v(t) ≤ vmax, x_con_vmax
    x(t0) == x0, initial_con
    m(tf) == mf, final_con
    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
    r(tf) → max

The (autonomous) optimal control problem is of the form:

    minimize  J(x, u, tf) = g(x(0), x(tf), tf)

    subject to

        ẋ(t) = f(x(t), u(t), tf), t in [0, tf] a.e.,

        ηl ≤ η(x(t), tf) ≤ ηu, 
        ϕl ≤ ϕ(x(0), x(tf), tf) ≤ ϕu, 

    where x(t) ∈ R³, u(t) ∈ R and tf ∈ R.

Declarations (* required):
╭────────┬────────┬──────────┬──────────┬───────────┬────────────┬─────────────╮
│ times* │ state* │ control* │ variable │ dynamics* │ objective* │ constraints │
├────────┼────────┼──────────┼──────────┼───────────┼────────────┼─────────────┤
│   V    │   V    │    V     │    V     │     V     │     V      │      V      │
╰────────┴────────┴──────────┴──────────┴───────────┴────────────┴─────────────╯

with

function F0(x)
    r, v, m = x
    D = Cd * v^2 * exp(-β*(r - 1))
    return [ v, -D/m - 1/r^2, 0 ]
end
function F1(x)
    r, v, m = x
    return [ 0, Tmax/m, -b*Tmax ]
end

You can plot the solution.

plot(prob.solution, size=(700, 700))

The solution is given by:

# initial costate
println("p0 = ", prob.solution.infos[:initial_costate])

# switching times
println("t1 = ", prob.solution.infos[:switching_times][1])
println("t2 = ", prob.solution.infos[:switching_times][2])
println("t3 = ", prob.solution.infos[:switching_times][3])

# final time
println("T  = ", prob.solution.infos[:final_time])

p0 = [3.945764658668555, 0.15039559623198723, 0.053712712939991955]
t1 = 0.023509684041475312
t2 = 0.059737380900899015
t3 = 0.10157134842460895
T  = 0.20204744057146434

References