Direct and indirect methods for the Goddard problem

Introduction

For this example, we consider the well-known Goddard problem^{[1] [2]} which 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_f$ 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

\[ r(t_f) \to \max\]

subject to the controlled 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.

We import the OptimalControl.jl package to define the optimal control problem and NLPModelsIpopt.jl to solve it. We import the Plots.jl package to plot the solution. The OrdinaryDiffEq.jl package is used to define the shooting function for the indirect method and the MINPACK.jl package permits to solve the shooting equation.

using OptimalControl  # to define the optimal control problem and more
using NLPModelsIpopt  # to solve the problem via a direct method
using OrdinaryDiffEq  # to get the Flow function from OptimalControl
using MINPACK         # NLE solver: use to solve the shooting equation
using Plots           # to plot the solution

Optimal control problem

We define the problem

const t0 = 0      # initial time
const r0 = 1      # initial altitude
const v0 = 0      # initial speed
const m0 = 1      # initial mass
const vmax = 0.1  # maximal authorized speed
const mf = 0.6    # final mass to target

goddard = @def begin # definition of the optimal control problem

    tf ∈ R, variable
    t ∈ [t0, tf], time
    x = (r, v, m) ∈ R³, state
    u ∈ R, control

    x(t0) == [r0, v0, m0]
    m(tf) == mf,         (1)
    0 ≤ u(t) ≤ 1
    r(t) ≥ r0
    0 ≤ v(t) ≤ vmax

    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))

    -r(tf) → min

end

# Dynamics
const Cd = 310
const Tmax = 3.5
const β = 500
const b = 2

F0(x) = begin
    r, v, m = x
    D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
    return [v, -D/m - 1/r^2, 0]
end

F1(x) = begin
    r, v, m = x
    return [0, Tmax/m, -b*Tmax]
end

Direct method

We then solve it

direct_sol = solve(goddard; grid_size=250)

▫ This is OptimalControl 2.0.1, solving with: collocation → adnlp → ipopt (cpu)

  📦 Configuration:
   ├─ Discretizer: collocation (grid_size = 250)
   ├─ Modeler: adnlp
   └─ Solver: ipopt

▫ This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:     4504
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     5259

Total number of variables............................:     1004
                     variables with only lower bounds:      251
                variables with lower and upper bounds:      501
                     variables with only upper bounds:        0
Total number of equality constraints.................:      754
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 -1.0100000e+00 9.00e-01 2.00e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -1.0091553e+00 8.99e-01 4.39e+01   1.1 1.19e+02    -  4.10e-03 8.29e-04f  1
   2 -1.0076086e+00 8.85e-01 3.41e+01   0.4 8.46e+00    -  2.02e-02 1.54e-02f  1
   3 -1.0015744e+00 7.58e-01 6.11e+02   0.5 6.85e+00    -  4.82e-01 1.44e-01f  1
   4 -1.0040637e+00 6.90e-01 1.28e+04   1.6 7.40e+00    -  6.88e-01 8.85e-02h  1
   5 -1.0054084e+00 6.10e-01 3.48e+04  -4.6 1.49e+00    -  1.82e-01 1.17e-01h  1
   6 -1.0056836e+00 5.93e-01 3.47e+04   1.7 7.88e+00    -  2.25e-01 2.72e-02f  1
   7 -1.0069013e+00 5.45e-01 1.39e+05   0.7 1.49e+00    -  2.86e-01 8.07e-02f  1
   8 -1.0092896e+00 4.27e-01 1.02e+05  -0.2 1.36e+00    -  2.79e-01 2.16e-01f  1
   9 -1.0120862e+00 6.66e-02 1.23e+05  -0.1 7.25e-01    -  3.54e-01 8.44e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -1.0082072e+00 1.19e-02 2.72e+04  -0.1 3.68e-01    -  6.09e-01 8.21e-01h  1
  11 -1.0072406e+00 1.32e-04 1.56e+04  -0.3 6.22e-01    -  7.62e-01 9.90e-01h  1
  12 -1.0071725e+00 1.27e-06 6.02e+02  -1.3 5.52e-02    -  9.75e-01 9.91e-01f  1
  13 -1.0071706e+00 1.30e-10 4.30e+02  -7.1 2.43e-04    -  9.89e-01 1.00e+00h  1
  14 -1.0071744e+00 2.57e-09 1.04e+03  -9.1 2.33e-03    -  9.59e-01 1.00e+00h  1
  15 -1.0072577e+00 1.17e-06 8.03e+02 -10.5 5.42e-02    -  2.29e-01 1.00e+00h  1
  16 -1.0077011e+00 4.13e-05 5.10e+02  -5.9 3.76e-01    -  2.60e-01 1.00e+00h  1
  17 -1.0078865e+00 4.11e-05 4.36e+02  -5.3 1.07e+00    -  3.14e-01 1.11e-01h  1
  18 -1.0083156e+00 2.01e-04 4.24e+03  -4.3 1.01e+01    -  1.00e+00 7.10e-02f  1
  19 -1.0088455e+00 1.37e-05 1.62e+03 -10.4 6.80e-02    -  6.57e-01 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -1.0106159e+00 5.71e-04 3.36e+02  -4.9 6.11e-01    -  1.00e+00 7.44e-01h  1
  21 -1.0099324e+00 1.58e-04 2.79e-03  -4.6 2.67e-01    -  1.00e+00 1.00e+00h  1
  22 -1.0107923e+00 9.18e-06 3.51e-03  -5.4 6.19e-02    -  1.00e+00 1.00e+00h  1
  23 -1.0122090e+00 3.87e-05 7.30e+01  -6.2 1.13e-01    -  1.00e+00 7.90e-01h  1
  24 -1.0124950e+00 5.99e-06 1.10e+00  -6.4 1.02e-01    -  1.00e+00 9.89e-01h  1
  25 -1.0125526e+00 3.20e-06 3.84e-05  -6.9 7.60e-02    -  1.00e+00 1.00e+00h  1
  26 -1.0125709e+00 1.92e-06 1.48e-05  -7.5 7.82e-02    -  1.00e+00 1.00e+00h  1
  27 -1.0125736e+00 6.35e-07 1.68e-01  -7.9 8.44e-02    -  1.00e+00 8.93e-01h  1
  28 -1.0125762e+00 3.57e-07 1.66e-06  -8.6 7.16e-02    -  1.00e+00 1.00e+00h  1
  29 -1.0125766e+00 2.93e-07 2.73e-07  -9.7 7.43e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30 -1.0125767e+00 5.86e-08 3.10e-08 -10.9 4.30e-02    -  1.00e+00 1.00e+00h  1
  31 -1.0125767e+00 1.86e-09 3.65e-10 -11.0 1.01e-02    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 31

                                   (scaled)                 (unscaled)
Objective...............:  -1.0125766794048943e+00   -1.0125766794048943e+00
Dual infeasibility......:   3.6466947942462423e-10    3.6466947942462423e-10
Constraint violation....:   1.0544083384189662e-09    1.8595086859196641e-09
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.1516659031614828e-11    2.1516659031614828e-11
Overall NLP error.......:   1.0544083384189662e-09    1.8595086859196641e-09


Number of objective function evaluations             = 32
Number of objective gradient evaluations             = 32
Number of equality constraint evaluations            = 32
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 32
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 31
Total seconds in IPOPT                               = 4.603

EXIT: Optimal Solution Found.

and plot the solution

plt = plot(direct_sol, label="direct", size=(800, 800))

Structure of the solution

We first determine visually the structure of the optimal solution which is composed of a bang arc with maximal control, followed by a singular arc, then by a boundary arc and the final arc is with zero control. In summary, the structure is: bang ($u=1$) → singular → boundary ($v=v_{\max}$) → off ($u=0$). Note that the switching function vanishes along the singular and boundary arcs.

t = time_grid(direct_sol)   # the time grid as a vector
x = state(direct_sol)       # the state as a function of time
u = control(direct_sol)     # the control as a function of time
p = costate(direct_sol)     # the costate as a function of time

H1 = Lift(F1)           # H1(x, p) = p' * F1(x)
φ(t) = H1(x(t), p(t))   # switching function
g(x) = vmax - x[2]      # state constraint v ≤ vmax

u_plot  = plot(t, u,     linewidth=2, label = "u(t)")
H1_plot = plot(t, φ,     linewidth=2, label = "H₁(x(t), p(t))")
g_plot  = plot(t, g ∘ x, linewidth=2, label = "g(x(t))")

plot(u_plot, H1_plot, g_plot, layout=(3,1), size=(600, 600))

We are now in position to solve the problem by an indirect shooting method. For an introduction to the indirect simple shooting method, see the Indirect simple shooting tutorial. We first define the four control laws in feedback form and their associated flows. For this we need to compute some Lie derivatives, namely Poisson brackets of Hamiltonians (themselves obtained as lifts to the cotangent bundle of vector fields), or derivatives of functions along a vector field. For instance, the control along the minimal order singular arcs is obtained as the quotient

\[u_s = -\frac{H_{001}}{H_{101}}\]

of length three Poisson brackets:

\[H_{001} = \{H_0,\{H_0,H_1\}\}, \quad H_{101} = \{H_1,\{H_0,H_1\}\}\]

where, for two Hamiltonians $H$ and $G$,

\[\{H,G\} := (\nabla_p H|\nabla_x G) - (\nabla_x H|\nabla_p G).\]

While the Lie derivative of a function $f$ wrt. a vector field $X$ is simply obtained as

\[(X \cdot f)(x) := f'(x) \cdot X(x),\]

and is used to the compute the control along the boundary arc,

\[u_b(x) = -(F_0 \cdot g)(x) / (F_1 \cdot g)(x),\]

as well as the associated multiplier for the order one state constraint on the velocity:

\[\mu(x, p) = H_{01}(x, p) / (F_1 \cdot g)(x).\]

With the help of differential geometry primitives, these expressions are straightforwardly translated into Julia code. For more details on the differential geometry tools, see Differential geometry tools. For a simpler example involving a minimal order singular arc, see Singular control.

# Controls
const u0 = 0                            # off control
const u1 = 1                            # bang control

H0 = Lift(F0)                           # H0(x, p) = p' * F0(x)
H01  = @Lie {H0, H1}
H001 = @Lie {H0, H01}
H101 = @Lie {H1, H01}
us(x, p) = -H001(x, p) / H101(x, p)     # singular control

ub(x) = -(F0⋅g)(x) / (F1⋅g)(x)          # boundary control
μ(x, p) = H01(x, p) / (F1⋅g)(x)         # multiplier associated to the state constraint g

# Flows
f0 = Flow(goddard, (x, p, tf) -> u0)
f1 = Flow(goddard, (x, p, tf) -> u1)
fs = Flow(goddard, (x, p, tf) -> us(x, p))
fb = Flow(goddard, (x, p, tf) -> ub(x), (x, u, tf) -> g(x), (x, p, tf) -> μ(x, p))

Shooting function

Then, we define the shooting function according to the optimal structure we have determined, that is a concatenation of four arcs. The shooting function has 7 unknowns (3 components of the initial costate p0 and 4 times: t1, t2, t3, tf) and 7 equations.

x0 = [r0, v0, m0] # initial state

function shoot!(s, p0, t1, t2, t3, tf)

    x1, p1 = f1(t0, x0, p0, t1)
    x2, p2 = fs(t1, x1, p1, t2)
    x3, p3 = fb(t2, x2, p2, t3)
    xf, pf = f0(t3, x3, p3, tf)

    s[1] = xf[3] - mf           # final mass constraint
    s[2:3] = pf[1:2] - [1, 0]   # transversality conditions: r(tf) and v(tf) are free
                                # objective is -r(tf) → p_r(tf) = 1, and v(tf) free → p_v(tf) = 0
    s[4] = H1(x1, p1)           # H1 = H01 = 0
    s[5] = H01(x1, p1)          # at the entrance of the singular arc
    s[6] = g(x2)                # g = 0 when entering the boundary arc
    s[7] = H0(xf, pf)           # since tf is free

end

Initial guess

To solve the problem by an indirect shooting method, we then need a good initial guess, that is a good approximation of the initial costate, the three switching times and the final time.

η = 1e-3
t13 = t[ abs.(φ.(t)) .≤ η ]
t23 = t[ 0 .≤ (g ∘ x).(t) .≤ η ]
p0 = p(t0)
t1 = min(t13...)
t2 = min(t23...)
t3 = max(t23...)
tf = t[end]

println("p0 = ", p0)
println("t1 = ", t1)
println("t2 = ", t2)
println("t3 = ", t3)
println("tf = ", tf)

# Norm of the shooting function at solution
using LinearAlgebra: norm
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")

p0 = [3.946831619547622, 0.14891030742499653, 0.05396551620088607]
t1 = 0.020197543560152617
t2 = 0.059784728938051736
t3 = 0.10179561954316918
tf = 0.20197543560152614

Norm of the shooting function: ‖s‖ = 1.8915305376704992

Indirect shooting

We aggregate the data to define the initial guess vector.

ξ_guess = [p0..., t1, t2, t3, tf] # initial guess

7-element Vector{Float64}:
 3.946831619547622
 0.14891030742499653
 0.05396551620088607
 0.020197543560152617
 0.059784728938051736
 0.10179561954316918
 0.20197543560152614

MINPACK.jl

We can use the MINPACK.jl package to solve the shooting equation. To compute the Jacobian of the shooting function we use the DifferentiationInterface.jl package with ForwardDiff.jl backend.

using DifferentiationInterface
import ForwardDiff
backend = AutoForwardDiff()

Let us define the problem to solve.

# auxiliary function with aggregated inputs
shoot!(s, ξ) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])

# Jacobian of the (auxiliary) shooting function
jshoot!(js, ξ) = jacobian!(shoot!, similar(ξ), js, backend, ξ)

We are now in position to solve the problem with the hybrj solver from MINPACK.jl through the fsolve function, providing the Jacobian. Let us solve the problem and retrieve the initial costate solution.

# resolution of S(ξ) = 0
minpack_sol = fsolve(shoot!, jshoot!, ξ_guess, show_trace=true)

# we retrieve the costate solution together with the times
ξ = minpack_sol.x
p0, t1, t2, t3, tf = ξ[1:3], ξ[4:end]...

println("MINPACK results:")
println("p0 = ", p0)
println("t1 = ", t1)
println("t2 = ", t2)
println("t3 = ", t3)
println("tf = ", tf)

# Norm of the shooting function at solution
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
@assert norm(s) < 1e-6 "Indirect shooting failed for MINPACK"

Iter     f(x) inf-norm    Step 2-norm      Step time
------   --------------   --------------   --------------
     1     1.889076e+00     0.000000e+00         0.284323
     2     6.350508e-02     3.864675e-03        53.922940
     3     2.263138e-02     1.051893e-03         0.069427
     4     8.512761e-03     7.573554e-04         0.053073
     5     1.089271e-03     1.621993e-07         0.034791
     6     1.731769e-03     6.487972e-07         0.033919
     7     5.837735e-05     1.459778e-06         0.034718
     8     7.127585e-05     6.305374e-08         0.034070
     9     2.930424e-05     2.262040e-09         0.034247
    10     5.289534e-07     8.272192e-11         0.043098
    11     6.800550e-07     2.274800e-12         0.035056
    12     6.250082e-08     4.748188e-14         0.034775
    13     1.958810e-09     2.216335e-15         0.035018
    14     1.649357e-09     3.440564e-17         0.034862
    15     2.821121e-10     2.134131e-19         0.034478
MINPACK results:
p0 = [3.9457646587813437, 0.15039559623311471, 0.053712712940263564]
t1 = 0.02350968404346425
t2 = 0.05973738089550592
t3 = 0.10157134842358748
tf = 0.20204744057152652

Norm of the shooting function: ‖s‖ = 2.8273316426846023e-10

NonlinearSolve.jl

Alternatively, we can use the NonlinearSolve.jl package to solve the shooting equation. The code is similar, but we use the solve function instead of fsolve. Let us define the problem.

shoot!(s, ξ, λ) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
prob = NonlinearProblem(shoot!, ξ_guess)

Now, let us solve the problem and retrieve the initial costate and times.

# resolution of S(ξ) = 0
sciml_sol = solve(prob; show_trace=Val(true))

# we retrieve the costate solution together with the times
ξ = sciml_sol.u
p0, t1, t2, t3, tf = ξ[1:3], ξ[4:end]...

println("\nNonlinearSolve results:")
println("p0 = ", p0)
println("t1 = ", t1)
println("t2 = ", t2)
println("t3 = ", t3)
println("tf = ", tf)

# Norm of the shooting function at solution
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
@assert norm(s) < 1e-6 "Indirect shooting failed for NonlinearSolve"

Algorithm: NewtonRaphson(
    descent = NewtonDescent(),
    autodiff = AutoForwardDiff(),
    vjp_autodiff = AutoReverseDiff(
        compile = false
    ),
    jvp_autodiff = AutoForwardDiff(),
    concrete_jac = Val{false}()
)

----    	-------------       	-----------
Iter    	f(u) inf-norm       	Step 2-norm
----    	-------------       	-----------
0       	1.88907616e+00      	0.00000000e+00
1       	6.35050789e-02      	6.21665070e-02
2       	5.18764339e-03      	6.17289944e-02
3       	4.50676595e-07      	6.84174864e-04
4       	1.62692082e-12      	6.54263450e-08
5       	3.31124017e-12      	3.33784233e-13
6       	1.12132525e-13      	8.52842174e-13
Final   	1.12132525e-13
----------------------

NonlinearSolve results:
p0 = [3.945764658689402, 0.15039559623165935, 0.05371271293970626]
t1 = 0.023509684041885595
t2 = 0.05973738089985927
t3 = 0.10157134842431804
tf = 0.20204744057100507

Norm of the shooting function: ‖s‖ = 1.2505769865390313e-13

Comparison of the two solvers

Let us compare the results obtained by the two solvers.

ξ_minpack = minpack_sol.x
ξ_sciml = sciml_sol.u

println("Comparison of the two solvers:")
println("MINPACK solution: ξ = ", ξ_minpack)
println("SciML solution:   ξ = ", ξ_sciml)
println("Difference:       Δξ = ", ξ_minpack - ξ_sciml)
println("Relative error:   ‖Δξ‖/‖ξ‖ = ", norm(ξ_minpack - ξ_sciml) / norm(ξ_minpack))

Comparison of the two solvers:
MINPACK solution: ξ = [3.9457646587813437, 0.15039559623311471, 0.053712712940263564, 0.02350968404346425, 0.05973738089550592, 0.10157134842358748, 0.20204744057152652]
SciML solution:   ξ = [3.945764658689402, 0.15039559623165935, 0.05371271293970626, 0.023509684041885595, 0.05973738089985927, 0.10157134842431804, 0.20204744057100507]
Difference:       Δξ = [9.194156547209786e-11, 1.455363607405502e-12, 5.57304202786213e-13, 1.5786538742901257e-12, -4.3533510130089326e-12, -7.305683835667764e-13, 5.214439990908204e-13]
Relative error:   ‖Δξ‖/‖ξ‖ = 2.327502608151321e-11

Benchmarking

The results found for by the two solvers are extremely close, so now, lets benchmark these two resolutions to compare their performances.

using BenchmarkTools

@benchmark fsolve(shoot!, jshoot!, ξ_guess; tol=1e-8, show_trace=false) #MINPACK

BenchmarkTools.Trial: 8 samples with 1 evaluation per sample.
 Range (min … max):  592.259 ms … 740.290 ms  ┊ GC (min … max):  7.02% … 22.95%
 Time  (median):     718.878 ms               ┊ GC (median):    22.26%
 Time  (mean ± σ):   688.240 ms ±  54.207 ms  ┊ GC (mean ± σ):  19.61% ±  6.07%

                                                      █          
  ▇▁▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▇▁▁▁▁▁▁▇ ▁
  592 ms           Histogram: frequency by time          740 ms <

 Memory estimate: 1.51 GiB, allocs estimate: 1646859.

@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve

BenchmarkTools.Trial: 3 samples with 1 evaluation per sample.
 Range (min … max):  2.066 s …    2.403 s  ┊ GC (min … max): 14.68% … 26.79%
 Time  (median):     2.372 s               ┊ GC (median):    25.44%
 Time  (mean ± σ):   2.280 s ± 186.410 ms  ┊ GC (mean ± σ):  22.66% ±  6.63%

  █                                                   █    █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁█ ▁
  2.07 s         Histogram: frequency by time          2.4 s <

 Memory estimate: 5.03 GiB, allocs estimate: 4469747.

According to the NonlinearSolve documentation, for small nonlinear systems, it could be faster to use the SimpleNewtonRaphson() descent algorithm.

@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false)) # NonlinearSolve

BenchmarkTools.Trial: 5 samples with 1 evaluation per sample.
 Range (min … max):  1.024 s …   1.245 s  ┊ GC (min … max):  6.96% … 23.97%
 Time  (median):     1.244 s              ┊ GC (median):    23.88%
 Time  (mean ± σ):   1.182 s ± 96.965 ms  ┊ GC (mean ± σ):  20.20% ±  7.38%

                                                          █  
  ▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  1.02 s         Histogram: frequency by time        1.25 s <

 Memory estimate: 2.87 GiB, allocs estimate: 1063930.

There exist different alternatives to solve the shooting equation, as shown in the benchmarks above.

Plot of the solution

We plot the solution of the indirect solution (in red) over the solution of the direct method (in blue).

f = f1 * (t1, fs) * (t2, fb) * (t3, f0) # concatenation of the flows
indirect_sol = f((t0, tf), x0, p0)      # compute the solution: state, costate, control...

plot!(plt, indirect_sol; label="indirect", color=2)

using Printf

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(direct_sol, indirect_sol)

    # get time grids
    t_dir = time_grid(direct_sol)
    t_ind = time_grid(indirect_sol)

    # create common time grid (union of both grids)
    t_common = unique(sort([t_dir..., t_ind...]))

    # interpolate both solutions on common grid
    x_dir_func = state(direct_sol)
    u_dir_func = control(direct_sol)
    x_ind_func = state(indirect_sol)
    u_ind_func = control(indirect_sol)

    x_dir = x_dir_func.(t_common)
    u_dir = u_dir_func.(t_common)
    x_ind = x_ind_func.(t_common)
    u_ind = u_ind_func.(t_common)

    v_dir = variable(direct_sol)
    o_dir = objective(direct_sol)
    i_dir = iterations(direct_sol)

    v_ind = variable(indirect_sol)
    o_ind = objective(indirect_sol)

    x_vars = ["r", "v", "m"]
    u_vars = ["u"]
    v_vars = ["tf"]

    println("┌─ Goddard problem: direct vs indirect")
    println("│")
    println("├─  Number of Iterations")
    @printf("│     Direct: %d\n", i_dir)
    println("│")

    # States
    println("├─  States (L2 Norms)")
    for i in eachindex(x_vars)
        xi_dir = [x_dir[k][i] for k in eachindex(t_common)]
        xi_ind = [x_ind[k][i] for k in eachindex(t_common)]
        L2_ae = L2_norm(t_common, xi_dir - xi_ind)
        L2_re = L2_ae / (0.5 * (L2_norm(t_common, xi_dir) + L2_norm(t_common, xi_ind)))
        @printf("│     %-6s Abs: %.3e   Rel: %.3e\n", x_vars[i], L2_ae, L2_re)
    end
    println("│")

    # Controls
    println("├─  Controls (L2 Norms)")
    for i in eachindex(u_vars)
        ui_dir = [u_dir[k][i] for k in eachindex(t_common)]
        ui_ind = [u_ind[k][i] for k in eachindex(t_common)]
        L2_ae = L2_norm(t_common, ui_dir - ui_ind)
        L2_re = L2_ae / (0.5 * (L2_norm(t_common, ui_dir) + L2_norm(t_common, ui_ind)))
        @printf("│     %-6s Abs: %.3e   Rel: %.3e\n", u_vars[i], L2_ae, L2_re)
    end
    println("│")

    # Variables
    println("├─  Variables")
    for i in eachindex(v_vars)
        vi_dir = v_dir[i]
        vi_ind = v_ind[i]
        vi_ae = abs(vi_dir - vi_ind)
        vi_re = vi_ae / (0.5 * (abs(vi_dir) + abs(vi_ind)))
        @printf("│     %-6s Abs: %.3e   Rel: %.3e\n", v_vars[i], vi_ae, vi_re)
    end
    println("│")

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

We now compare numerically the direct and indirect solutions. The comparison includes L2 norms for the state and control variables, as well as absolute and relative errors for the variable (final time) and the objective. The L2 norm is computed using the trapezoidal rule over the time grid of the direct solution.

print_numerical_comparisons(direct_sol, indirect_sol)

┌─ Goddard problem: direct vs indirect
│
├─  Number of Iterations
│     Direct: 31
│
├─  States (L2 Norms)
│     r      Abs: 3.810e-07   Rel: 8.418e-07
│     v      Abs: 1.030e-05   Rel: 3.281e-04
│     m      Abs: 1.790e-05   Rel: 5.852e-05
│
├─  Controls (L2 Norms)
│     u      Abs: 4.817e-03   Rel: 2.373e-02
│
├─  Variables
│     tf     Abs: 7.200e-05   Rel: 3.564e-04
│
├─  Objective
│            Abs: 3.548e-07   Rel: 3.504e-07
└─

Hamiltonian verification

The Hamiltonian should be constant along the optimal trajectory and equal to zero since the final time is free. Let us verify this for both the direct and indirect solutions.

# Hamiltonian function
H(x, p, u) = H0(x, p) + u * H1(x, p)

# Direct solution
t_dir = time_grid(direct_sol)
x_dir = state(direct_sol)
u_dir = control(direct_sol)
p_dir = costate(direct_sol)
H_direct = [H(x_dir(ti), p_dir(ti), u_dir(ti)) for ti in t_dir]

# Indirect solution
t_ind = time_grid(indirect_sol)
x_ind = state(indirect_sol)
u_ind = control(indirect_sol)
p_ind = costate(indirect_sol)
H_indirect = [H(x_ind(ti), p_ind(ti), u_ind(ti)) for ti in t_ind]

println("Direct method:")
println("  H(t0) = ", H_direct[1])
println("  H variation: max|H(t) - H(t0)| = ", maximum(abs.(H_direct .- H_direct[1])))

println("\nIndirect method:")
println("  H(t0) = ", H_indirect[1])
println("  H variation: max|H(t) - H(t0)| = ", maximum(abs.(H_indirect .- H_indirect[1])))

Direct method:
  H(t0) = -0.005482855660857239
  H variation: max|H(t) - H(t0)| = 0.005964115097448206

Indirect method:
  H(t0) = 1.204508714991448e-12
  H variation: max|H(t) - H(t0)| = 3.304895982061766e-10

We can also plot the Hamiltonian along the trajectory to verify its constancy visually.

# Plot
plot(t_dir, H_direct, label="Direct", linewidth=2)
plot!(t_ind, H_indirect, label="Indirect", linewidth=2, linestyle=:dash)
xlabel!("Time")
ylabel!("H(t)")
title!("Hamiltonian along the trajectory")

References