Robert Hutchings Goddard was an American engineer, professor, physicist, and inventor who is credited with creating and building the world's first liquid-fueled rocket which he successfully launched on March 16, 1926. Goddard and his team launched 34 rockets between 1926 and 1941, achieving altitudes as high as 2.6 km (1.6 mi) and speeds as high as 885 km/h (550 mph).
From Wikipedia
This problem was drawn from the COPS3 benchmark.
Before we even talk about optimization, we will create a simple model of a rocket in flight: state variables, controls, and dynamics.
We will take a very simplified view and consider everything to be one-dimensional: our rocket will go straight up (and straight down). So our first two state variables that vary across time are
Our rocket is burning fuel as it goes, getting lighter and lighter in the process. Therefore we need a third state variable
We will start at altitude $h_0$, velocity 0, and mass $m_0$. We will define $m_c$ to be fraction of rocket that is the rocket itself. $m(t)$ is thus $\geq m_c m_0$.
We will have one control for our rocket: the thrust $T(t)$. We will assume that it can be changed rapidly, which is not completely unreasonable for a liquid-fuel rocket. There is however a maximum thrust of $T_{max}$.
We need to consider what forces act on our rocket other than thrust. The most obvious perhaps is gravity
$$g(h) = g_0 \left( \frac{h(0)}{h} \right)^2$$which is a function of height. The other key force is drag
$$D(h,v) = D_c v^2 exp\left( -h_c \left( \frac{h-h(0)}{h(0)} \right) \right)$$which is a complicated function of velocity (squared), height (to adjust for air density), and a drag coefficient $D_c$ determined by the rocket shape.
Given all the above, we can link together everything with three equations to control the dynamics of the rocket:
where $c$ is a constant measuring the impulse of the rocket fuel.
The goal is maximize the final altitude of our rocket. The only thing we can change, other than the rocket properties itself, are the way we vary our thrust over time. We will formulate an optimization problem to find the optimal thrust profile.
We will use a discretized model of time, with a fixed number of time steps, $n$. We will make the time step size $\Delta t$, and thus the final time $t_f = n \cdot \Delta t$, a variable in the problem. We can now write out objective as $$ \max h(t_f) $$
Given this, we now need to discretize our dynamic equations. For example, how do we relate $v(t)$ and $v(t+\Delta t)$? To do so we will use numerical integration. One simple idea is the "rectangular rule", e.g. $$ f(t+\Delta t) - f(t) = \Delta t v(t) $$
We will use the more accurate but more complicated trapezoidal rule (other options are possible!). In particular, in the context our problem we have that $$ f(t+\Delta t) - f(t) = \frac{ v(t + \Delta t) + v(t) }{ 0.5 \Delta t } $$
The three equations that relate the states, controls, and forces correspond to constraints in the optimization problem.
For the purposes of this example we will use dimension free (scaled) parameters. In particular:
Finally we will use $n$ time steps.
We will be using an interior point solver. You can think of it as basically a gradient descent inside the space of feasible solutions. We can help out the solver by providing a good initial solution that is not wildly different from the final solution. We'll use the following solution:
# Load JuMP, and a nonlinear solver Ipopt
# Other options we have for this problem class
# is KNITRO.jl, a commercial solver.
using JuMP, Ipopt
# Create JuMP model, using Ipopt as the solver
#mod = Model(solver=IpoptSolver(print_level=0))
mod = Model(solver=IpoptSolver())
# Constants
# Note that all parameters in the model have been normalized
# to be dimensionless. See the COPS3 paper for more info.
h_0 = 1 # Initial height
v_0 = 0 # Initial velocity
m_0 = 1 # Initial mass
g_0 = 1 # Gravity at the surface
# Parameters
T_c = 3.5 # Used for thrust
h_c = 500 # Used for drag
v_c = 620 # Used for drag
m_c = 0.6 # Fraction of initial mass left at end
# Derived parameters
c = 0.5*sqrt(g_0*h_0) # Thrust-to-fuel mass
m_f = m_c*m_0 # Final mass
D_c = 0.5*v_c*m_0/g_0 # Drag scaling
T_max = T_c*g_0*m_0 # Maximum thrust
# Time steps
n = 800
800
# Set up variables and constraints
# Time step with initial guess
@defVar(mod, Δt ≥ 0, start = 1/n)
# Store a useful subexpression, "time of flight"
@defNLExpr(t_f, Δt*n)
# State variables
@defVar(mod, v[0:n] ≥ 0) # Velocity
@defVar(mod, h[0:n] ≥ h_0) # Height
@defVar(mod, m_f ≤ m[0:n] ≤ m_0) # Mass
# Control: thrust
@defVar(mod, 0 ≤ T[0:n] ≤ T_max)
# Provide starting solution
# Could have done this at same time as @defVar
for k in 0:n
setValue(h[k], 1)
setValue(v[k], (k/n)*(1 - (k/n)))
setValue(m[k], (m_f - m_0)*(k/n) + m_0)
setValue(T[k], T_max/2)
end
# Objective: maximize altitude at end of time of flight
@setObjective(mod, Max, h[n])
# Initial conditions
@addConstraint(mod, v[0] == v_0)
@addConstraint(mod, h[0] == h_0)
@addConstraint(mod, m[0] == m_0)
@addConstraint(mod, m[n] == m_f)
# Forces
# We'll define these as expressions too
# Wherever they appear, they will effectively be
# replaced by these longer expressions, keeps them
# nice and clean
# Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )
@defNLExpr(drag[j=0:n], D_c*(v[j]^2)*exp(-h_c*(h[j]-h_0)/h_0))
# Grav(h) = g0 * (h0 / h)^2
@defNLExpr(grav[j=0:n], g_0*(h_0/h[j])^2)
# Dynamics
for j in 1:n
# h' = v
# Rectangular integration
# @addNLConstraint(mod, h[j] == h[j-1] + Δt*v[j-1])
# Trapezoidal integration
@addNLConstraint(mod,
h[j] == h[j-1] + 0.5*Δt*(v[j]+v[j-1]))
# v' = (T-D(h,v))/m - g(h)
# Rectangular integration
# @addNLConstraint(mod, v[j] == v[j-1] + Δt*(
# (T[j-1] - drag[j-1])/m[j-1] - grav[j-1]))
# Trapezoidal integration
@addNLConstraint(mod,
v[j] == v[j-1] + 0.5*Δt*(
(T[j ] - drag[j ] - m[j ]*grav[j ])/m[j ] +
(T[j-1] - drag[j-1] - m[j-1]*grav[j-1])/m[j-1] ))
# m' = -T/c
# Rectangular integration
# @addNLConstraint(mod, m[j] == m[j-1] - Δt*T[j-1]/c)
# Trapezoidal integration
@addNLConstraint(mod,
m[j] == m[j-1] - 0.5*Δt*(T[j] + T[j-1])/c)
end
# Solve for the control and state
println("Solving...")
status = solve(mod)
# Display results
println("Solver status: ", status)
println("Max height: ", getObjectiveValue(mod))
Solving... This is Ipopt version 3.12.1, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 15204 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 48800 Total number of variables............................: 3205 variables with only lower bounds: 1603 variables with lower and upper bounds: 1602 variables with only upper bounds: 0 Total number of equality constraints.................: 2404 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 3.50e-02 1.83e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 -1.0476357e+00 3.10e-03 7.60e+03 -1.0 2.38e-01 2.0 9.18e-01 9.67e-01f 1 2 -1.0330218e+00 1.03e-03 8.61e+04 -1.0 2.31e+00 - 3.71e-01 7.53e-01h 1 3 -1.0108102e+00 3.89e-04 1.79e+05 -1.0 1.27e+00 - 1.00e+00 6.40e-01h 1 4 -1.0058561e+00 1.58e-04 2.97e+05 -1.0 6.84e-01 - 1.00e+00 9.69e-01h 1 5 -1.0078548e+00 9.05e-05 2.16e+05 -1.0 8.08e-01 - 1.00e+00 6.33e-01f 1 6 -1.0079653e+00 1.75e-05 4.11e+04 -1.0 4.89e-01 - 1.00e+00 1.00e+00f 1 7 -1.0078191e+00 1.38e-06 1.02e+02 -1.0 1.70e-01 - 1.00e+00 1.00e+00f 1 8 -1.0078115e+00 3.32e-10 2.15e+02 -1.7 3.11e-03 - 1.00e+00 1.00e+00h 1 9 -1.0078157e+00 2.53e-10 4.46e-01 -1.7 1.20e-02 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 -1.0078162e+00 6.07e-12 1.15e+01 -3.8 4.17e-04 - 1.00e+00 1.00e+00h 1 11 -1.0078877e+00 1.06e-07 5.24e-03 -3.8 5.52e-02 - 1.00e+00 1.00e+00f 1 12 -1.0078881e+00 2.27e-13 1.65e-06 -3.8 3.09e-04 - 1.00e+00 1.00e+00h 1 13 -1.0079601e+00 1.06e-07 2.95e+02 -5.7 5.49e-02 - 9.80e-01 1.00e+00f 1 14 -1.0095526e+00 8.48e-05 1.22e+00 -5.7 2.23e+00 - 1.00e+00 9.65e-01h 1 15 -1.0112897e+00 5.88e-05 2.95e+00 -5.7 1.78e+00 - 1.00e+00 9.90e-01f 1 16 -1.0111071e+00 2.22e-06 2.07e-02 -5.7 7.35e-01 - 1.00e+00 1.00e+00f 1 17 -1.0111190e+00 6.83e-09 6.56e-05 -5.7 5.17e-02 - 1.00e+00 1.00e+00h 1 18 -1.0111191e+00 4.12e-14 1.20e-09 -5.7 3.56e-04 - 1.00e+00 1.00e+00h 1 19 -1.0122683e+00 2.30e-05 2.99e+01 -8.6 8.98e-01 - 7.15e-01 8.77e-01f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 20 -1.0127029e+00 1.51e-05 8.14e+00 -8.6 8.55e-01 - 7.50e-01 7.80e-01h 1 21 -1.0128032e+00 1.01e-05 3.13e+00 -8.6 1.26e+00 - 6.59e-01 7.31e-01h 1 22 -1.0128269e+00 5.32e-06 1.22e+00 -8.6 1.40e+00 - 6.46e-01 7.33e-01h 1 23 -1.0128326e+00 2.58e-06 3.96e-01 -8.6 1.38e+00 - 6.98e-01 7.71e-01f 1 24 -1.0128339e+00 1.13e-06 1.87e-02 -8.6 1.24e+00 - 9.40e-01 8.79e-01f 1 25 -1.0128341e+00 2.77e-07 5.12e-05 -8.6 9.73e-01 - 1.00e+00 1.00e+00f 1 26 -1.0128341e+00 3.97e-09 1.13e-06 -8.6 2.34e-01 - 1.00e+00 1.00e+00h 1 27 -1.0128341e+00 1.03e-10 7.14e-09 -8.6 3.13e-02 - 1.00e+00 1.00e+00h 1 Number of Iterations....: 27 (scaled) (unscaled) Objective...............: -1.0128340615485458e+00 -1.0128340615485458e+00 Dual infeasibility......: 7.1350593378980852e-09 7.1350593378980852e-09 Constraint violation....: 1.0326103860869296e-10 1.0326103860869296e-10 Complementarity.........: 2.5081820874400327e-09 2.5081820874400327e-09 Overall NLP error.......: 7.1350593378980852e-09 7.1350593378980852e-09 Number of objective function evaluations = 28 Number of objective gradient evaluations = 28 Number of equality constraint evaluations = 28 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 28 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 27 Total CPU secs in IPOPT (w/o function evaluations) = 0.259 Total CPU secs in NLP function evaluations = 0.265 EXIT: Optimal Solution Found. Solver status: Optimal Max height: 1.0128340615485458
# Can visualize the state and control variables
using Gadfly
h_plot = plot(x=(0:n)*getValue(Δt),y=getValue(h)[1:n], Geom.line,
Guide.xlabel("Time (s)"), Guide.ylabel("Altitude"))
m_plot = plot(x=(0:n)*getValue(Δt),y=getValue(m)[1:n], Geom.line,
Guide.xlabel("Time (s)"), Guide.ylabel("Mass"))
v_plot = plot(x=(0:n)*getValue(Δt),y=getValue(v)[1:n], Geom.line,
Guide.xlabel("Time (s)"), Guide.ylabel("Velocity"))
T_plot = plot(x=(0:n)*getValue(Δt),y=getValue(T)[1:n], Geom.line,
Guide.xlabel("Time (s)"), Guide.ylabel("Thrust"))
draw(SVG(6inch, 6inch), vstack( hstack(h_plot,m_plot),
hstack(v_plot,T_plot)))
WARNING: groupby(xs, keyfunc) should be groupby(keyfunc, xs)