Rocket Launch
We will solve an optimal control problem that aims to minimize final time.
Problem Statement and Model
In this problem, we seek to minimize the final time $t_f$ of a rocket launch between 0 and 10 metres. Minimizing fuel use is also important for this problem. The force- or thrust- of the rocket can be between -1.1 and 1.1, while the velocity of the rocket cannot exceed 1.7. The drag resistance of the rocket is proportional to the square of the velocity. The mass of the rocket decreases as fuel is burned over time. Our optimization guidelines are subject to these constraints:
\[\begin{aligned} &&\min t_f \\ &&&\frac{ds}{dt}= v \\ &&&\frac{dv}{dt}= (u-0.2*(v^2))/2\\ &&&\frac{dm}{dt}= -0.01(u^2)\\ &&&0 \leq v(t) \leq 1.7 \\ &&&-1.1 \leq u(t) \leq 1.1 \\ &&&s(0) = 0 \\ &&&v(0) = 0 \\ &&&m(0) = 1 \\ &&&s(t_f) = 10 \\ &&&v(t_f) = 0 \\ \end{aligned}\]
where $v(t)$ represents the rocket's velocity, $s(t)$ represents the rocket's position and $m(t)$ represents the rocket's mass. Additionally, $u(t)$ represents the rocket's force (or thrust).
Modeling in InfiniteOpt
First we must import $InfiniteOpt$ and other packages.
using JuMP, InfiniteOpt, Ipopt, Plots;Model Initialization
We then initialize the infinite model with InfiniteModel, selecting Ipopt as the desired optimizer to solve this model.
model = InfiniteModel(Ipopt.Optimizer);Infinite Parameter Definition
We are defining the infinite parameter $t \in [0, 1]$. The bounds give the optimizer parameters to iterate over.
@infinite_parameter(model, t in [0,1], num_supports = 50, derivative_method = FiniteDifference(Backward()));Infinite Variable Definition
Moving on to specifying variables, we have v as velocity of the rocket, u as the force of the rocket. The position, p, and mass, m, vary over time. The final time, t_f, must be between 0.1 and 100 seconds.
@variable(model, 0 <= v <= 1.7, Infinite(t))
@variable(model, -1.1 <= u <= 1.1, Infinite(t))
@variable(model, 0.1 <= t_f <= 100)
@variable(model, p, Infinite(t))
@variable(model, m, Infinite(t));Initial and Final Condition Definition
Now that the model is defined, we then want to specify both our initial and final conditions.
@constraint(model, v(0) == 0)
@constraint(model, v(1) == 0)
@constraint(model, p(0) == 0)
@constraint(model, p(1) == 10)
@constraint(model, m(0) == 1)
@constraint(model, m(1) >= 0.2);Objective Definition
We now add the objective, using @objective, which is to minimize final time t_f.
@objective(model, Min, t_f);Constraint Definition
We finally are adding our constraints, which are defining the ODEs for our system model.
@constraint(model, v * t_f == deriv(p, t))
@constraint(model, ((u - (0.2 * (v ^ 2))) * t_f) == (deriv(v, t)) * m)
@constraint(model, (-0.01 * (u ^ 2)) * t_f == deriv(m, t));Problem Solution
Now that everything is set up, we can solve the model by using @optimize!
optimize!(model)This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.
Number of nonzeros in equality constraint Jacobian...: 996
Number of nonzeros in inequality constraint Jacobian.: 1
Number of nonzeros in Lagrangian Hessian.............: 600
Total number of variables............................: 351
variables with only lower bounds: 0
variables with lower and upper bounds: 101
variables with only upper bounds: 0
Total number of equality constraints.................: 302
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 1
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.0999999e-01 1.00e+01 1.58e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.0009999e-01 1.00e+01 1.43e+02 -1.0 2.90e+03 -4.0 1.41e-02 9.73e-05H 1
2 1.1725017e-01 9.84e+00 5.08e+06 -1.0 5.01e+03 - 7.75e-05 1.64e-02f 1
3 1.1789798e-01 9.83e+00 5.08e+06 -1.0 1.85e+03 - 5.50e-04 2.16e-04h 2
4 1.1798656e-01 9.83e+00 5.08e+06 -1.0 2.61e+03 - 1.73e-03 2.93e-05h 4
5 1.1805455e-01 9.83e+00 5.08e+06 -1.0 5.17e+03 - 9.29e-04 2.06e-05h 4
6 1.1867019e-01 9.83e+00 5.08e+06 -1.0 6.32e+03 - 1.54e-04 1.54e-04s 10
7r 1.1867019e-01 9.83e+00 9.99e+02 1.0 0.00e+00 - 0.00e+00 0.00e+00R 1
8r 1.1791305e-01 9.83e+00 9.94e+02 1.0 1.82e+02 - 7.58e-03 5.06e-03f 1
9r 1.1550385e-01 9.82e+00 9.80e+02 1.0 3.38e+01 - 6.09e-02 1.45e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r 1.0853227e-01 9.71e+00 8.62e+02 1.0 1.16e+01 - 4.52e-01 1.20e-01f 1
11r 1.3689437e-01 8.80e+00 1.80e+03 1.0 1.55e+01 - 1.80e-01 8.00e-01f 1
12 1.4610978e-01 8.79e+00 8.28e+02 -1.0 3.47e+04 - 3.66e-04 1.01e-03f 2
13 1.4890163e-01 8.79e+00 7.94e+02 -1.0 3.04e+04 - 1.23e-03 2.91e-04h 3
14 1.4960075e-01 8.79e+00 7.93e+02 -1.0 2.09e+04 - 2.09e-03 7.93e-05h 5
15 1.4982760e-01 8.79e+00 7.93e+02 -1.0 1.41e+04 - 4.33e-03 2.75e-05h 7
16 1.5013079e-01 8.79e+00 7.93e+02 -1.0 5.79e+03 - 6.44e-03 3.95e-05h 7
17 1.6404713e-01 8.77e+00 3.19e+03 -1.0 9.86e+03 - 1.21e-03 1.72e-03f 2
18 1.7461253e-01 8.76e+00 2.51e+03 -1.0 1.12e+04 - 1.64e-03 1.24e-03f 2
19 1.7775177e-01 8.76e+00 2.49e+03 -1.0 7.46e+03 - 3.64e-03 3.82e-04h 3
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.7846508e-01 8.76e+00 2.49e+03 -1.0 5.43e+03 - 6.26e-03 8.76e-05h 5
21 1.7976322e-01 8.76e+00 2.49e+03 -1.0 3.14e+03 - 4.55e-03 1.58e-04h 5
22 2.1717764e-01 8.72e+00 2.09e+03 -1.0 7.85e+03 - 1.58e-03 4.24e-03w 1
23 2.1899940e-01 8.72e+00 2.09e+03 -1.0 1.29e+03 - 6.57e-03 2.27e-04w 1
24 2.3996486e-01 8.70e+00 2.07e+03 -1.0 3.24e+03 - 2.72e-03 2.55e-03w 1
25 1.9847043e-01 8.74e+00 2.18e+03 -1.0 2.39e+03 - 1.58e-03 2.12e-03f 1
26 2.0370158e-01 8.73e+00 2.17e+03 -1.0 4.69e+03 - 5.57e-03 6.24e-04h 3
27 2.0868653e-01 8.73e+00 2.17e+03 -1.0 3.28e+03 - 8.37e-03 6.02e-04h 3
28 2.2593378e-01 8.71e+00 2.17e+03 -1.0 2.53e+03 - 6.48e-03 2.08e-03h 2
29 2.4744655e-01 8.69e+00 2.15e+03 -1.0 3.05e+03 - 3.63e-03 2.55e-03h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 2.7725933e-01 8.66e+00 2.13e+03 -1.0 2.14e+03 - 1.09e-02 3.61e-03h 1
31 3.0489726e-01 8.63e+00 2.13e+03 -1.0 1.34e+03 - 1.27e-02 3.37e-03h 1
32 3.4019435e-01 8.59e+00 2.11e+03 -1.0 1.14e+03 - 7.86e-03 4.29e-03h 2
33 3.8627668e-01 8.54e+00 2.10e+03 -1.0 8.31e+02 - 2.45e-02 5.64e-03h 2
34 4.3593059e-01 8.49e+00 2.09e+03 -1.0 3.47e+02 - 3.54e-02 6.12e-03h 2
35 5.3565226e-01 8.38e+00 2.06e+03 -1.0 4.27e+02 - 2.73e-02 1.23e-02h 2
36 6.5465550e-01 8.26e+00 2.03e+03 -1.0 2.27e+02 - 7.53e-02 1.50e-02h 2
37 8.9104567e-01 8.01e+00 1.96e+03 -1.0 9.16e+01 - 7.45e-02 3.02e-02h 2
38 1.2985846e+00 7.58e+00 1.86e+03 -1.0 6.99e+01 - 2.23e-01 5.40e-02h 2
39 2.1415554e+00 6.67e+00 1.64e+03 -1.0 2.04e+01 - 4.79e-01 1.19e-01h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 4.2026938e+00 4.42e+00 1.09e+03 -1.0 9.30e+00 - 9.42e-01 3.37e-01h 2
41 8.0245946e+00 3.24e-01 1.19e+01 -1.0 6.25e+00 -1.8 1.00e+00 9.90e-01h 1
42 7.9987268e+00 2.13e-02 7.02e+01 -1.0 2.52e+00 - 9.30e-01 9.90e-01h 1
43 8.6845877e+00 2.90e-01 2.83e+04 -1.0 3.47e+00 - 7.57e-01 1.00e+00f 1
44 9.9285930e+00 5.21e-02 4.54e+05 -1.0 3.53e+00 - 9.54e-01 1.00e+00H 1
45 1.0034156e+01 5.41e-03 1.05e-03 -1.0 5.18e-01 - 1.00e+00 1.00e+00h 1
46 8.6094459e+00 3.48e-01 1.46e+06 -3.8 1.88e+00 - 8.53e-01 1.00e+00f 1
47 7.6023969e+00 4.98e-01 3.75e+05 -3.8 4.31e+00 - 7.43e-01 1.00e+00h 1
48 7.4901387e+00 1.17e-01 8.80e+04 -3.8 2.78e+00 - 7.66e-01 8.01e-01h 1
49 7.4674801e+00 3.92e-02 1.57e+04 -3.8 1.05e+00 - 8.22e-01 6.68e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 7.4556190e+00 1.50e-03 2.26e+03 -3.8 1.02e+00 - 8.56e-01 1.00e+00h 1
51 7.4554953e+00 2.76e-05 2.65e-05 -3.8 1.41e-01 - 1.00e+00 1.00e+00h 1
52 7.4488628e+00 7.23e-05 9.30e+02 -5.7 1.06e-01 - 9.85e-01 9.22e-01h 1
53 7.4483375e+00 6.04e-07 2.30e-07 -5.7 8.51e-03 - 1.00e+00 1.00e+00h 1
54 7.4482497e+00 1.32e-08 3.65e-01 -8.6 1.07e-03 - 1.00e+00 9.98e-01h 1
55 7.4482495e+00 5.68e-14 7.77e-14 -8.6 2.50e-06 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 55
(scaled) (unscaled)
Objective...............: 7.4482495165303710e+00 7.4482495165303710e+00
Dual infeasibility......: 7.7656631125577746e-14 7.7656631125577746e-14
Constraint violation....: 5.6843418860808015e-14 5.6843418860808015e-14
Variable bound violation: 3.2404592624652052e-39 3.2404592624652052e-39
Complementarity.........: 2.5059601286904100e-09 2.5059601286904100e-09
Overall NLP error.......: 2.5059601286904100e-09 2.5059601286904100e-09
Number of objective function evaluations = 148
Number of objective gradient evaluations = 53
Number of equality constraint evaluations = 148
Number of inequality constraint evaluations = 148
Number of equality constraint Jacobian evaluations = 57
Number of inequality constraint Jacobian evaluations = 57
Number of Lagrangian Hessian evaluations = 55
Total seconds in IPOPT = 0.043
EXIT: Optimal Solution Found.
Extract and Plot the Results
We can now extract the optimized results and plot them to visualize four results graphically: Position over time, velocity over time, mass over time and force over time. In addition, we can print the final time value of t_f.
tf_data = value(t_f)
t_data = supports(t) .* tf_data
v_data = value(v)
p_data = value(p)
m_data = value(m)
u_data = value(u);Printing the final time value.
println("Minimized Final Time: ", tf_data," seconds");Minimized Final Time: 7.448249516530371 seconds
Creating the four plots as described above.
plot1 = plot(t_data, v_data, title="Velocity Over Time", linecolor = :red, xlabel="Time", label = "Velocity", ylabel="Velocity", lw=2);
plot2 = plot(t_data, p_data, title="Position Over Time", linecolor = :orange, xlabel="Time", label = "Position", ylabel="Position", lw=2);
plot3 = plot(t_data, m_data, title="Mass Over Time", linecolor = :yellow, xlabel="Time", label = "Position", ylabel="Position", lw=2);
plot4 = plot(t_data, u_data, title="Force Over Time", linecolor = :limegreen, xlabel="Time", label = "Position", ylabel="Position", lw=2);Visualizing all four plots on one display.
display(plot(plot1, plot2, plot3, plot4, layout=(4,1), size = (1000, 1000)))Maintenance Tests
These are here to ensure this example stays up to date.
using Test
tol = 1E-4
@test termination_status(model) == MOI.LOCALLY_SOLVED
@test has_values(model)
@test isapprox(objective_value(model), 7.4482495165303710, atol=tol)Test PassedThis page was generated using Literate.jl.