{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Disclaimer\n",
"This notebook is only working under the versions:\n",
"\n",
"- JuMP 0.19 (unreleased, but currently in master)\n",
"\n",
"- MathOptInterface 0.4.1\n",
"\n",
"- Ipopt 0.4.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Description**: Shows how to solve a nonlinear rocketry control problem.\n",
"\n",
"**Author**: Iain Dunning\n",
"\n",
"**License**: ![\"Creative](\"https://i.creativecommons.org/l/by-sa/4.0/88x31.png\")
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Rocket Control with JuMP\n",
"\n",
"*This problem was drawn from the [COPS3](http://www.mcs.anl.gov/~more/cops/cops3.pdf) benchmark.*\n",
"\n",
"Our goal is to maximize the final altitude of a vertically launched rocket. We can control the thrust of the rocket, and must take account of the rocket mass, fuel consumption rate, gravity, and aerodynamic drag.\n",
"\n",
"Let us consider the basic description of the model (for the full description, including parameters for the rocket, see the COPS3 PDF)\n",
"\n",
"### Overview\n",
"\n",
"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. To approximate the derivatives in the problem we will use the [trapezoidal rule](http://en.wikipedia.org/wiki/Trapezoidal_rule).\n",
"\n",
"### State and Control\n",
"\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$. Our goal is thus to maximize $h(t_f)$. Each of these corresponds to a JuMP variable indexed by the time step.\n",
"\n",
"### Dynamics\n",
"\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. 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, and for convenience we will represent the forces with nonlinear expressions that we define seperately with `@NLexpression`."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuMP, Ipopt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"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.0128340648308067\n"
]
}
],
"source": [
"# Create JuMP model, using Ipopt as the solver\n",
"mod = Model(optimizer=IpoptOptimizer(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",
"# Parameters\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",
"# Derived parameters\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(mod, 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(mod, Max, h[n])\n",
"\n",
"# Initial conditions\n",
"@constraints(mod, 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(mod, 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(mod, grav[j=1:n], g_0*(h_0/h[j])^2)\n",
"# Time of flight\n",
"@NLexpression(mod, t_f, Δt*n)\n",
"\n",
"# Dynamics\n",
"for j in 2:n\n",
" # h' = v\n",
" # Rectangular integration\n",
" # @NLconstraint(mod, h[j] == h[j-1] + Δt*v[j-1])\n",
" # Trapezoidal integration\n",
" @NLconstraint(mod,\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",
" # Rectangular integration\n",
" # @NLconstraint(mod, v[j] == v[j-1] + Δt*(\n",
" # (T[j-1] - drag[j-1])/m[j-1] - grav[j-1]))\n",
" # Trapezoidal integration\n",
" @NLconstraint(mod,\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",
" # Rectangular integration\n",
" # @NLconstraint(mod, m[j] == m[j-1] - Δt*T[j-1]/c)\n",
" # Trapezoidal integration\n",
" @NLconstraint(mod,\n",
" m[j] == m[j-1] - 0.5*Δt*(T[j] + T[j-1])/c)\n",
"end\n",
"\n",
"# Provide starting solution\n",
"# for k in 0:n\n",
"# setvalue(h[k], 1)\n",
"# setvalue(v[k], (k/n)*(1 - (k/n)))\n",
"# setvalue(m[k], (m_f - m_0)*(k/n) + m_0)\n",
"# setvalue(T[k], T_max/2)\n",
"# end\n",
"\n",
"# Solve for the control and state\n",
"println(\"Solving...\")\n",
"status = JuMP.optimize(mod)\n",
"\n",
"# Display results\n",
"# println(\"Solver status: \", status)\n",
"println(\"Max height: \", JuMP.objectivevalue(mod))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Can visualize the state and control variables\n",
"using Gadfly"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Compose.SVG(152.39999999999998mm, 152.39999999999998mm, IOBuffer(data=UInt8[...], readable=true, writable=true, seekable=true, append=false, size=59365, maxsize=Inf, ptr=59366, mark=-1), nothing, \"img-bf9a5313\", 0, Compose.SVGPropertyFrame[], Dict{Type,Union{Compose.Property, Void}}(), DataStructures.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.485mm, 81.2mm), (147.4mm, 81.2mm), (147.4mm, 133.115mm), (92.485mm, 133.115mm)])=>\"img-bf9a5313-4\",Compose.ClipPrimitive{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}(Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}[(20.305mm, 81.2mm), (71.2mm, 81.2mm), (71.2mm, 133.115mm), (20.305mm, 133.115mm)])=>\"img-bf9a5313-17\",Compose.ClipPrimitive{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}(Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}[(94.8317mm, 5.0mm), (147.4mm, 5.0mm), (147.4mm, 56.915mm), (94.8317mm, 56.915mm)])=>\"img-bf9a5313-30\",Compose.ClipPrimitive{Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}}(Tuple{Measures.Length{:mm,Float64},Measures.Length{:mm,Float64}}[(21.9783mm, 5.0mm), (71.2mm, 5.0mm), (71.2mm, 56.915mm), (21.9783mm, 56.915mm)])=>\"img-bf9a5313-43\"), Tuple{Compose.FormPrimitive,String}[], Set{AbstractString}(), true, false, nothing, true, \"img-bf9a5313-52\", false, 52, Set(AbstractString[\"C:\\\\Users\\\\joaquimgarcia\\\\.julia\\\\v0.6\\\\Gadfly\\\\src\\\\gadfly.js\"]), Set(Tuple{AbstractString,AbstractString}[(\"Snap.svg\", \"Snap\"), (\"Gadfly\", \"Gadfly\")]), AbstractString[\"fig.select(\\\"#img-bf9a5313-5\\\")\\n .init_gadfly();\", \"fig.select(\\\"#img-bf9a5313-18\\\")\\n .init_gadfly();\", \"fig.select(\\\"#img-bf9a5313-31\\\")\\n .init_gadfly();\", \"fig.select(\\\"#img-bf9a5313-44\\\")\\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.008845317510043124, 1.0156261087342169, 0.2176906350200863, -0.016252217468433822, 0.0mm, 0.0mm, 0.0mm, 0.0mm), Compose.IdentityTransform()))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h_plot = plot(x=(1:n)*JuMP.resultvalue.(Δt),y=JuMP.resultvalue.(h)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Altitude\"))\n",
"m_plot = plot(x=(1:n)*JuMP.resultvalue.(Δt),y=JuMP.resultvalue.(m)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Mass\"))\n",
"v_plot = plot(x=(1:n)*JuMP.resultvalue.(Δt),y=JuMP.resultvalue.(v)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Velocity\"))\n",
"T_plot = plot(x=(1:n)*JuMP.resultvalue.(Δt),y=JuMP.resultvalue.(T)[:], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Thrust\"))\n",
"draw(SVG(6inch, 6inch), vstack( hstack(h_plot,m_plot),\n",
" hstack(v_plot,T_plot)))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}