{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"title: Rocket Control\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Originally Contributed by**: Iain Dunning"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This tutorial shows how to solve a nonlinear rocketry control problem.\n",
"The problem was drawn from the [COPS3](http://www.mcs.anl.gov/~more/cops/cops3.pdf) benchmark."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our goal is to maximize the final altitude of a vertically launched rocket. \n",
"We can control the thrust of the rocket, and must take account of \n",
"the rocket mass, fuel consumption rate, gravity, and aerodynamic drag."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us consider the basic description of the model (for the full description, \n",
"including parameters for the rocket, see the COPS3 PDF)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Overview \n",
"We will use a discretized model of time, with a fixed number of time steps, $n$. \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. \n",
"To approximate the derivatives in the problem we will use the [trapezoidal rule](http://en.wikipedia.org/wiki/Trapezoidal_rule)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### State and Control \n",
"We will have three state variables:\n",
"\n",
"* Velocity, $v$\n",
"* Altitude, $h$\n",
"* Mass of rocket and remaining fuel, $m$\n",
"\n",
"and a single control variable, thrust $T$. \n",
"Our goal is thus to maximize $h(t_f)$. \n",
"Each of these corresponds to a JuMP variable indexed by the time step."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Dynamics\n",
"We have three equations that control the dynamics of the rocket:\n",
"\n",
"Rate of ascent: $$h^\\prime = v$$\n",
"Acceleration: $$v^\\prime = \\frac{T - D(h,v)}{m} - g(h)$$\n",
"Rate of mass loss: $$m^\\prime = -\\frac{T}{c}$$\n",
"\n",
"where drag $D(h,v)$ is a function of altitude and velocity, and gravity $g(h)$ is a function of altitude."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These forces are defined as\n",
"\n",
"$$D(h,v) = D_c v^2 exp\\left( -h_c \\left( \\frac{h-h(0)}{h(0)} \\right) \\right)$$\n",
"and\n",
"$$g(h) = g_0 \\left( \\frac{h(0)}{h} \\right)^2$$\n",
"\n",
"The three rate equations correspond to JuMP constraints, \n",
"and for convenience we will represent the forces with nonlinear expressions."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving...\n",
"\n",
"******************************************************************************\n",
"This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
" Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
" For more information visit http://projects.coin-or.org/Ipopt\n",
"******************************************************************************\n",
"\n",
"Max height: 1.0128340648308058\n"
]
}
],
"source": [
"using JuMP, Ipopt\n",
"\n",
"# Create JuMP model, using Ipopt as the solver\n",
"rocket = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0))\n",
"\n",
"# Constants\n",
"# Note that all parameters in the model have been normalized\n",
"# to be dimensionless. See the COPS3 paper for more info.\n",
"h_0 = 1 # Initial height\n",
"v_0 = 0 # Initial velocity\n",
"m_0 = 1 # Initial mass\n",
"g_0 = 1 # Gravity at the surface\n",
"\n",
"T_c = 3.5 # Used for thrust\n",
"h_c = 500 # Used for drag\n",
"v_c = 620 # Used for drag\n",
"m_c = 0.6 # Fraction of initial mass left at end\n",
"\n",
"c = 0.5 * sqrt(g_0 * h_0) # Thrust-to-fuel mass\n",
"m_f = m_c * m_0 # Final mass\n",
"D_c = 0.5 * v_c * m_0 / g_0 # Drag scaling\n",
"T_max = T_c * g_0 * m_0 # Maximum thrust\n",
"\n",
"n = 800 # Time steps\n",
"\n",
"@variables(rocket, begin\n",
" Δt ≥ 0, (start = 1/n) # Time step\n",
" # State variables\n",
" v[1:n] ≥ 0 # Velocity\n",
" h[1:n] ≥ h_0 # Height\n",
" m_f ≤ m[1:n] ≤ m_0 # Mass\n",
" # Control\n",
" 0 ≤ T[1:n] ≤ T_max # Thrust\n",
"end)\n",
"\n",
"# Objective: maximize altitude at end of time of flight\n",
"@objective(rocket, Max, h[n])\n",
"\n",
"# Initial conditions\n",
"@constraints(rocket, begin\n",
" v[1] == v_0\n",
" h[1] == h_0\n",
" m[1] == m_0\n",
" m[n] == m_f\n",
"end)\n",
"\n",
"# Forces\n",
"# Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )\n",
"@NLexpression(rocket, drag[j = 1:n], D_c * (v[j]^2) * exp(-h_c * (h[j] - h_0) / h_0))\n",
"# Grav(h) = go * (h0 / h)^2\n",
"@NLexpression(rocket, grav[j = 1:n], g_0 * (h_0 / h[j])^2)\n",
"# Time of flight\n",
"@NLexpression(rocket, t_f, Δt * n)\n",
"\n",
"# Dynamics\n",
"for j in 2:n\n",
" # h' = v\n",
" \n",
" # Rectangular integration\n",
" # @NLconstraint(rocket, h[j] == h[j - 1] + Δt * v[j - 1])\n",
" \n",
" # Trapezoidal integration\n",
" @NLconstraint(rocket,\n",
" h[j] == h[j - 1] + 0.5 * Δt * (v[j] + v[j - 1]))\n",
"\n",
" # v' = (T-D(h,v))/m - g(h)\n",
" \n",
" # Rectangular integration\n",
" # @NLconstraint(rocket, v[j] == v[j - 1] + Δt *(\n",
" # (T[j - 1] - drag[j - 1]) / m[j - 1] - grav[j - 1]))\n",
" \n",
" # Trapezoidal integration\n",
" @NLconstraint(rocket,\n",
" v[j] == v[j-1] + 0.5 * Δt * (\n",
" (T[j] - drag[j] - m[j] * grav[j]) / m[j] +\n",
" (T[j - 1] - drag[j - 1] - m[j - 1] * grav[j - 1]) / m[j - 1]))\n",
"\n",
" # m' = -T/c\n",
"\n",
" # Rectangular integration\n",
" # @NLconstraint(rocket, m[j] == m[j - 1] - Δt * T[j - 1] / c)\n",
" \n",
" # Trapezoidal integration\n",
" @NLconstraint(rocket,\n",
" m[j] == m[j - 1] - 0.5 * Δt * (T[j] + T[j-1]) / c)\n",
"end\n",
"\n",
"# Solve for the control and state\n",
"println(\"Solving...\")\n",
"status = optimize!(rocket)\n",
"\n",
"# Display results\n",
"# println(\"Solver status: \", status)\n",
"println(\"Max height: \", objective_value(rocket))"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Can visualize the state and control variables\n",
"using Gadfly"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"SVG(152.39999999999998mm, 152.39999999999998mm, IOBuffer(data=UInt8[...], readable=true, writable=true, seekable=true, append=false, size=61717, maxsize=Inf, ptr=61718, mark=-1), nothing, \"img-a1f9bcce\", 0, Compose.SVGPropertyFrame[], Dict{Type,Union{Nothing, Compose.Property}}(Compose.Property{Compose.JSIncludePrimitive} => nothing,Compose.Property{Compose.FontSizePrimitive} => nothing,Compose.Property{Compose.StrokePrimitive} => nothing,Compose.Property{Compose.FontPrimitive} => nothing,Compose.Property{Compose.SVGAttributePrimitive} => nothing,Compose.Property{Compose.SVGClassPrimitive} => nothing,Compose.Property{Compose.FillPrimitive} => nothing,Compose.Property{Compose.LineWidthPrimitive} => nothing,Compose.Property{Compose.JSCallPrimitive} => nothing,Compose.Property{Compose.StrokeDashPrimitive} => nothing…), OrderedCollections.OrderedDict{Compose.ClipPrimitive,String}(Compose.ClipPrimitive{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}(Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}[(92.48499999999999mm, 81.19999999999999mm), (147.39999999999998mm, 81.19999999999999mm), (147.39999999999998mm, 133.11499999999998mm), (92.48499999999999mm, 133.11499999999998mm)]) => \"img-a1f9bcce-4\",Compose.ClipPrimitive{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}(Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}[(20.305mm, 81.19999999999999mm), (71.19999999999999mm, 81.19999999999999mm), (71.19999999999999mm, 133.11499999999998mm), (20.305mm, 133.11499999999998mm)]) => \"img-a1f9bcce-41\",Compose.ClipPrimitive{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}(Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}[(94.83166666666665mm, 5.0mm), (147.39999999999998mm, 5.0mm), (147.39999999999998mm, 56.91499999999999mm), (94.83166666666665mm, 56.91499999999999mm)]) => \"img-a1f9bcce-76\",Compose.ClipPrimitive{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}(Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}[(21.97833333333333mm, 5.0mm), (71.19999999999999mm, 5.0mm), (71.19999999999999mm, 56.91499999999999mm), (21.97833333333333mm, 56.91499999999999mm)]) => \"img-a1f9bcce-115\"), Tuple{Compose.FormPrimitive,String}[], Set{AbstractString}(), true, false, nothing, true, \"img-a1f9bcce-140\", false, 140, Set(AbstractString[\"/home/mbesancon/.julia/packages/Gadfly/USbaq/src/gadfly.js\"]), Set(Tuple{AbstractString,AbstractString}[(\"Snap.svg\", \"Snap\"), (\"Gadfly\", \"Gadfly\")]), AbstractString[\"fig.select(\\\"#img-a1f9bcce-5\\\")\\n .init_gadfly();\", \"fig.select(\\\"#img-a1f9bcce-42\\\")\\n .init_gadfly();\", \"fig.select(\\\"#img-a1f9bcce-77\\\")\\n .init_gadfly();\", \"fig.select(\\\"#img-a1f9bcce-116\\\")\\n .init_gadfly();\"], false, :none, (Measures.BoundingBox{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}},Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}((21.97833333333333mm, 5.0mm), (49.22166666666666mm, 51.91499999999999mm)), Compose.UnitBox{Float64,Float64,Float64,Float64}(-0.008845317510043122, 1.0156261087342169, 0.21769063502008623, -0.016252217468433822, 0.0mm, 0.0mm, 0.0mm, 0.0mm), Compose.IdentityTransform()))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h_plot = plot(x = (1:n) * value.(Δt), y = value.(h)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Altitude\"))\n",
"m_plot = plot(x = (1:n) * value.(Δt), y = value.(m)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Mass\"))\n",
"v_plot = plot(x = (1:n) * value.(Δt), y = value.(v)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Velocity\"))\n",
"T_plot = plot(x = (1:n) * value.(Δt), y = value.(T)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Thrust\"))\n",
"draw(SVG(6inch, 6inch), vstack(hstack(h_plot, m_plot), hstack(v_plot, T_plot)))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.1",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}