Goddard Rocket Control

using Gridap, PDENLPModels, PDEOptimizationProblems

This tutorial shows how to solve a nonlinear rocketry control problem. The problem was drawn from the COPS3 benchmark, and was used in JuMP's tutorials rocket control.

PDEOptimizationProblems contains a list of predefined PDE-constrained optimization problems, e.g. Goddard rocket control problem via the function rocket. Note that the models in PDEOptimizationProblems are in Lagrangian form, so the implementation considers a Mayer to Lagrange formulation that adds an additional (constant) function H := h(T) see in the code.

T = 0.2 # final time
n = 800 # number of cells

nlp = rocket(n = n, T = T)

GridapPDENLPModel
  Problem name: Goddard Rocket n=800
   All variables: ████████████████████ 4800   All constraints: ████████████████████ 3200  
            free: ████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 801               free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
           lower: ███████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 1600             lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
           upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
         low/upp: ██████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 2399           low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
           fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                fixed: ████████████████████ 3200  
          infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
            nnzh: ( 99.24% sparsity)   87971           linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
                                                    nonlinear: ████████████████████ 3200  
                                                         nnzj: ( 99.58% sparsity)   64001 

  Counters:
             obj: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 grad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 cons: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
        cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0             cons_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
           jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                  jac: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0              jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
         jac_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0            jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
       jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jtprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0           jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
      jtprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                 hess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0                hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     
           jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0               jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0     

The return model is a GridapPDENLPModel and, in particular, is an instance of an AbstractNLPModel that can be solved with any tool from the JuliaSmoothOptimizers organization. Then, using NLPModelsIpopt.jl, we solve the model and get the required information.

using NLPModelsIpopt

stats = ipopt(nlp, x0 = nlp.meta.x0)
obj_ipopt, con_ipopt = stats.objective, stats.primal_feas
(hh, Hh, vh, mh), uh = split_vectors(nlp, stats.solution)

(Base.Generator{UnitRange{Int64}, PDENLPModels.var"#26#32"{Vector{Float64}, Gridap.MultiField.MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{Float64}}, PDENLPModels.var"#sum_old#30", PDENLPModels.var"#leng#28"}}(PDENLPModels.var"#26#32"{Vector{Float64}, Gridap.MultiField.MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{Float64}}, PDENLPModels.var"#sum_old#30", PDENLPModels.var"#leng#28"}([1.0000000782063823, 1.0000003130096415, 1.0000007047709742, 1.0000012538374403, 1.0000019605418038, 1.0000028252023891, 1.0000038481229483, 1.000005029592539, 1.0000063698854187, 1.0000078692609484  …  0.00029721802875868556, 0.00029721802875868556, 0.0003279979227687752, 0.0003279979227687752, 0.0003276886489833572, 0.0003276886489833572, 0.0004048727303939332, 0.0004048727303939332, 0.0004043900185087497, 0.0004043900185087497], MultiFieldFESpace(), PDENLPModels.var"#sum_old#30"(), PDENLPModels.var"#leng#28"()), 1:4), [3.499593124492011, 3.499593124492011, 3.499583233072741, 3.499583233072741, 3.4995730612885647, 3.4995730612885647, 3.499562551364897, 3.499562551364897, 3.4995516768785424, 3.4995516768785424  …  0.00029721802875868556, 0.00029721802875868556, 0.0003279979227687752, 0.0003279979227687752, 0.0003276886489833572, 0.0003276886489833572, 0.0004048727303939332, 0.0004048727303939332, 0.0004043900185087497, 0.0004043900185087497])

Finally, we can plot the functions, and the results match JuMP's tutorial and COPS 3 report.

using Plots
gr()

h₀, m₀, mᵪ = 1.0, 1.0, 0.6
p = Plots.plot(
  Plots.plot(0:T/n:T, vcat(h₀, hh); xlabel = "Time (s)", ylabel = "Altitude"),
  Plots.plot(0:T/n:T, vcat(m₀, mh, mᵪ * m₀); xlabel = "Time (s)", ylabel = "Mass"),
  Plots.plot(0:T/n:T, vcat(0., vh); xlabel = "Time (s)", ylabel = "Velocity"),
  Plots.plot(0:T/(length(uh) - 1):T, uh; xlabel = "Time (s)", ylabel = "Thrust"),
  layout = (2, 2),
  legend = false,
  margin = 1Plots.cm,
)

An alternative is also to intepolate the entire solution over the domain as an FEFunction (Gridap's function type) and save the interpolation in a VTK file.

Ypde = nlp.pdemeta.Ypde
h_h  = FEFunction(Ypde.spaces[1], hh)
H_h  = FEFunction(Ypde.spaces[2], Hh)
v_h  = FEFunction(Ypde.spaces[3], vh)
m_h  = FEFunction(Ypde.spaces[4], mh)
u_h  = FEFunction(nlp.pdemeta.Ycon, uh)

writevtk(
  nlp.pdemeta.tnrj.trian,
  "results_robot",
  cellfields=["uh" => u_h, "hh" => h_h, "Hh" => H_h, "vh" => v_h, "mh" => m_h],
)

(["results_robot.vtu"],)