{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Rocket Control with JuMP\n",
"\n",
"\n",
"\n",
"> 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).\n",
"\n",
"*From [Wikipedia](https://en.wikipedia.org/wiki/Robert_H._Goddard)*\n",
"\n",
"*This problem was drawn from the [COPS3](http://www.mcs.anl.gov/~more/cops/cops3.pdf) benchmark.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simple Rocket Flight Model\n",
"\n",
"Before we even talk about optimization, we will create a simple model of a rocket in flight: state variables, controls, and dynamics.\n",
"\n",
"### State variables\n",
"\n",
"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 \n",
"\n",
"* **velocity**, $v(t)$\n",
"* **altitude**, $h(t)$\n",
"\n",
"Our rocket is burning fuel as it goes, getting lighter and lighter in the process. Therefore we need a third state variable\n",
"\n",
"* **Mass of rocket and remaining fuel**, $m(t)$\n",
"\n",
"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$.\n",
"\n",
"### Control\n",
"\n",
"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}$.\n",
"\n",
"### Dynamics\n",
"\n",
"We need to consider what forces act on our rocket other than thrust. The most obvious perhaps is **gravity**\n",
"\n",
"$$g(h) = g_0 \\left( \\frac{h(0)}{h} \\right)^2$$\n",
"\n",
"which is a function of height. The other key force is **drag**\n",
"\n",
"$$D(h,v) = D_c v^2 exp\\left( -h_c \\left( \\frac{h-h(0)}{h(0)} \\right) \\right)$$\n",
"\n",
"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.\n",
"\n",
"Given all the above, we can link together everything with three equations to 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 $c$ is a constant measuring the impulse of the rocket fuel."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optimization problem\n",
"\n",
"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.\n",
"\n",
"### Dealing with time\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. We can now write out objective as\n",
"$$\n",
"\\max h(t_f)\n",
"$$\n",
"\n",
"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.\n",
"$$\n",
"f(t+\\Delta t) - f(t) = \\Delta t v(t)\n",
"$$\n",
"\n",
"We will use the more accurate but more complicated [trapezoidal rule](http://en.wikipedia.org/wiki/Trapezoidal_rule) (other options are possible!). In particular, in the context our problem we have that\n",
"$$\n",
"f(t+\\Delta t) - f(t) = \\frac{ v(t + \\Delta t) + v(t) }{ 0.5 \\Delta t }\n",
"$$\n",
"\n",
"The three equations that relate the states, controls, and forces correspond to **constraints** in the optimization problem.\n",
"\n",
"### Data\n",
"\n",
"For the purposes of this example we will use dimension free (scaled) parameters. In particular:\n",
"\n",
"* $T_{max} = 3.5 g_0 m_0$\n",
"* $D_c = 0.5 v_c \\frac{m_0}{g_0}$\n",
"* $c = 0.5 (g_0 h_0 )^{0.5}$\n",
"* So $h_0 = m_0 = g_0 = 1$\n",
"* and we choose $h_c = 500$, $m_c = 0.6$, $v_c = 620$\n",
"\n",
"Finally we will use $n$ time steps.\n",
"\n",
"### Initial Solution\n",
"\n",
"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:\n",
"\n",
"* $t_f = n \\Delta t = 1$\n",
"* $h(t) = 1$\n",
"* $v(t) = \\frac{t}{t_f}(1 - \\frac{t}{t_f})$\n",
"* $m(t) = (m_f - m_0)\\frac{t}{t_f} + m_0$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Building the model"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Load JuMP, and a nonlinear solver Ipopt\n",
"# Other options we have for this problem class\n",
"# is KNITRO.jl, a commercial solver.\n",
"using JuMP, Ipopt"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"$$ \\begin{alignat*}{1}\\min\\quad & 0\\\\\n",
"\\text{Subject to} \\quad\\end{alignat*}\n",
" $$"
],
"text/plain": [
"Feasibility problem with:\n",
" * 0 linear constraints\n",
" * 0 variables\n",
"Solver set to Ipopt"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Create JuMP model, using Ipopt as the solver\n",
"#mod = Model(solver=IpoptSolver(print_level=0))\n",
"mod = Model(solver=IpoptSolver())"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"800"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 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",
"# Time steps\n",
"n = 800"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Set up variables and constraints\n",
"\n",
"# Time step with initial guess\n",
"@defVar(mod, Δt ≥ 0, start = 1/n)\n",
"# Store a useful subexpression, \"time of flight\"\n",
"@defNLExpr(t_f, Δt*n)\n",
"\n",
"# State variables\n",
"@defVar(mod, v[0:n] ≥ 0) # Velocity\n",
"@defVar(mod, h[0:n] ≥ h_0) # Height\n",
"@defVar(mod, m_f ≤ m[0:n] ≤ m_0) # Mass\n",
"\n",
"# Control: thrust\n",
"@defVar(mod, 0 ≤ T[0:n] ≤ T_max)\n",
"\n",
"# Provide starting solution\n",
"# Could have done this at same time as @defVar\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",
"# Objective: maximize altitude at end of time of flight\n",
"@setObjective(mod, Max, h[n])\n",
"\n",
"# Initial conditions\n",
"@addConstraint(mod, v[0] == v_0)\n",
"@addConstraint(mod, h[0] == h_0)\n",
"@addConstraint(mod, m[0] == m_0)\n",
"@addConstraint(mod, m[n] == m_f)\n",
"\n",
"# Forces\n",
"# We'll define these as expressions too\n",
"# Wherever they appear, they will effectively be\n",
"# replaced by these longer expressions, keeps them\n",
"# nice and clean\n",
"\n",
"# Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )\n",
"@defNLExpr(drag[j=0:n], D_c*(v[j]^2)*exp(-h_c*(h[j]-h_0)/h_0))\n",
"# Grav(h) = g0 * (h0 / h)^2\n",
"@defNLExpr(grav[j=0:n], g_0*(h_0/h[j])^2)\n",
"\n",
"# Dynamics\n",
"for j in 1:n\n",
" # h' = v\n",
" # Rectangular integration\n",
" # @addNLConstraint(mod, h[j] == h[j-1] + Δt*v[j-1])\n",
" # Trapezoidal integration\n",
" @addNLConstraint(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",
" # @addNLConstraint(mod, v[j] == v[j-1] + Δt*(\n",
" # (T[j-1] - drag[j-1])/m[j-1] - grav[j-1]))\n",
" # Trapezoidal integration\n",
" @addNLConstraint(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",
" # @addNLConstraint(mod, m[j] == m[j-1] - Δt*T[j-1]/c)\n",
" # Trapezoidal integration\n",
" @addNLConstraint(mod,\n",
" m[j] == m[j-1] - 0.5*Δt*(T[j] + T[j-1])/c)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving...\n",
"This is Ipopt version 3.12.1, running with linear solver mumps.\n",
"NOTE: Other linear solvers might be more efficient (see Ipopt documentation).\n",
"\n",
"Number of nonzeros in equality constraint Jacobian...: 15204\n",
"Number of nonzeros in inequality constraint Jacobian.: 0\n",
"Number of nonzeros in Lagrangian Hessian.............: 48800\n",
"\n",
"Total number of variables............................: 3205\n",
" variables with only lower bounds: 1603\n",
" variables with lower and upper bounds: 1602\n",
" variables with only upper bounds: 0\n",
"Total number of equality constraints.................: 2404\n",
"Total number of inequality constraints...............: 0\n",
" inequality constraints with only lower bounds: 0\n",
" inequality constraints with lower and upper bounds: 0\n",
" inequality constraints with only upper bounds: 0\n",
"\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 0 -1.0100000e+00 3.50e-02 1.83e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n",
" 1 -1.0476357e+00 3.10e-03 7.60e+03 -1.0 2.38e-01 2.0 9.18e-01 9.67e-01f 1\n",
" 2 -1.0330218e+00 1.03e-03 8.61e+04 -1.0 2.31e+00 - 3.71e-01 7.53e-01h 1\n",
" 3 -1.0108102e+00 3.89e-04 1.79e+05 -1.0 1.27e+00 - 1.00e+00 6.40e-01h 1\n",
" 4 -1.0058561e+00 1.58e-04 2.97e+05 -1.0 6.84e-01 - 1.00e+00 9.69e-01h 1\n",
" 5 -1.0078548e+00 9.05e-05 2.16e+05 -1.0 8.08e-01 - 1.00e+00 6.33e-01f 1\n",
" 6 -1.0079653e+00 1.75e-05 4.11e+04 -1.0 4.89e-01 - 1.00e+00 1.00e+00f 1\n",
" 7 -1.0078191e+00 1.38e-06 1.02e+02 -1.0 1.70e-01 - 1.00e+00 1.00e+00f 1\n",
" 8 -1.0078115e+00 3.32e-10 2.15e+02 -1.7 3.11e-03 - 1.00e+00 1.00e+00h 1\n",
" 9 -1.0078157e+00 2.53e-10 4.46e-01 -1.7 1.20e-02 - 1.00e+00 1.00e+00h 1\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 10 -1.0078162e+00 6.07e-12 1.15e+01 -3.8 4.17e-04 - 1.00e+00 1.00e+00h 1\n",
" 11 -1.0078877e+00 1.06e-07 5.24e-03 -3.8 5.52e-02 - 1.00e+00 1.00e+00f 1\n",
" 12 -1.0078881e+00 2.27e-13 1.65e-06 -3.8 3.09e-04 - 1.00e+00 1.00e+00h 1\n",
" 13 -1.0079601e+00 1.06e-07 2.95e+02 -5.7 5.49e-02 - 9.80e-01 1.00e+00f 1\n",
" 14 -1.0095526e+00 8.48e-05 1.22e+00 -5.7 2.23e+00 - 1.00e+00 9.65e-01h 1\n",
" 15 -1.0112897e+00 5.88e-05 2.95e+00 -5.7 1.78e+00 - 1.00e+00 9.90e-01f 1\n",
" 16 -1.0111071e+00 2.22e-06 2.07e-02 -5.7 7.35e-01 - 1.00e+00 1.00e+00f 1\n",
" 17 -1.0111190e+00 6.83e-09 6.56e-05 -5.7 5.17e-02 - 1.00e+00 1.00e+00h 1\n",
" 18 -1.0111191e+00 4.12e-14 1.20e-09 -5.7 3.56e-04 - 1.00e+00 1.00e+00h 1\n",
" 19 -1.0122683e+00 2.30e-05 2.99e+01 -8.6 8.98e-01 - 7.15e-01 8.77e-01f 1\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 20 -1.0127029e+00 1.51e-05 8.14e+00 -8.6 8.55e-01 - 7.50e-01 7.80e-01h 1\n",
" 21 -1.0128032e+00 1.01e-05 3.13e+00 -8.6 1.26e+00 - 6.59e-01 7.31e-01h 1\n",
" 22 -1.0128269e+00 5.32e-06 1.22e+00 -8.6 1.40e+00 - 6.46e-01 7.33e-01h 1\n",
" 23 -1.0128326e+00 2.58e-06 3.96e-01 -8.6 1.38e+00 - 6.98e-01 7.71e-01f 1\n",
" 24 -1.0128339e+00 1.13e-06 1.87e-02 -8.6 1.24e+00 - 9.40e-01 8.79e-01f 1\n",
" 25 -1.0128341e+00 2.77e-07 5.12e-05 -8.6 9.73e-01 - 1.00e+00 1.00e+00f 1\n",
" 26 -1.0128341e+00 3.97e-09 1.13e-06 -8.6 2.34e-01 - 1.00e+00 1.00e+00h 1\n",
" 27 -1.0128341e+00 1.03e-10 7.14e-09 -8.6 3.13e-02 - 1.00e+00 1.00e+00h 1\n",
"\n",
"Number of Iterations....: 27\n",
"\n",
" (scaled) (unscaled)\n",
"Objective...............: -1.0128340615485458e+00 -1.0128340615485458e+00\n",
"Dual infeasibility......: 7.1350593378980852e-09 7.1350593378980852e-09\n",
"Constraint violation....: 1.0326103860869296e-10 1.0326103860869296e-10\n",
"Complementarity.........: 2.5081820874400327e-09 2.5081820874400327e-09\n",
"Overall NLP error.......: 7.1350593378980852e-09 7.1350593378980852e-09\n",
"\n",
"\n",
"Number of objective function evaluations = 28\n",
"Number of objective gradient evaluations = 28\n",
"Number of equality constraint evaluations = 28\n",
"Number of inequality constraint evaluations = 0\n",
"Number of equality constraint Jacobian evaluations = 28\n",
"Number of inequality constraint Jacobian evaluations = 0\n",
"Number of Lagrangian Hessian evaluations = 27\n",
"Total CPU secs in IPOPT (w/o function evaluations) = 0.259\n",
"Total CPU secs in NLP function evaluations = 0.265\n",
"\n",
"EXIT: Optimal Solution Found.\n",
"Solver status: Optimal\n",
"Max height: 1.0128340615485458\n"
]
}
],
"source": [
"# Solve for the control and state\n",
"println(\"Solving...\")\n",
"status = solve(mod)\n",
"\n",
"# Display results\n",
"println(\"Solver status: \", status)\n",
"println(\"Max height: \", getObjectiveValue(mod))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Can visualize the state and control variables\n",
"using Gadfly"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"SVG(152.39999999999998,152.39999999999998,IOBuffer(data=Uint8[...], readable=true, writable=true, seekable=true, append=false, size=54579, maxsize=Inf, ptr=54580, mark=-1),nothing,\"fig-2a34be1329b044df9921e3cf270e22e3\",0,[],Dict{Type{T<:Top},Union(Nothing,Property{P<:PropertyPrimitive})}(),[ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(21.97833333333333,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(21.97833333333333,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-42\",ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(20.305,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(20.305,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-18\",ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(94.83166666666665,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(94.83166666666665,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-30\",ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(92.48499999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(92.48499999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-6\"],Set{String}({}),true,false,nothing,true,\"fig-2a34be1329b044df9921e3cf270e22e3-element-48\",false,48,String[\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\",\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\",\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\",\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\"],(String,String)[(\"Snap.svg\",\"Snap\"),(\"Gadfly\",\"Gadfly\"),(\"Gadfly\",\"Gadfly\"),(\"Gadfly\",\"Gadfly\"),(\"Gadfly\",\"Gadfly\")],String[\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-5\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\",\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-17\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\",\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-29\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\",\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-41\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\"],false,:none)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: groupby(xs, keyfunc) should be groupby(keyfunc, xs)\n"
]
}
],
"source": [
"h_plot = plot(x=(0:n)*getValue(Δt),y=getValue(h)[1:n], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Altitude\"))\n",
"m_plot = plot(x=(0:n)*getValue(Δt),y=getValue(m)[1:n], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Mass\"))\n",
"v_plot = plot(x=(0:n)*getValue(Δt),y=getValue(v)[1:n], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Velocity\"))\n",
"T_plot = plot(x=(0:n)*getValue(Δt),y=getValue(T)[1:n], 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)))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.3.7",
"language": "julia",
"name": "julia 0.3"
},
"language_info": {
"name": "julia",
"version": "0.3.7"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
![](\"https://upload.wikimedia.org/wikipedia/commons/7/7c/Goddard_and_Rocket.jpg\")
\n",
"\n",
"> 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).\n",
"\n",
"*From [Wikipedia](https://en.wikipedia.org/wiki/Robert_H._Goddard)*\n",
"\n",
"*This problem was drawn from the [COPS3](http://www.mcs.anl.gov/~more/cops/cops3.pdf) benchmark.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simple Rocket Flight Model\n",
"\n",
"Before we even talk about optimization, we will create a simple model of a rocket in flight: state variables, controls, and dynamics.\n",
"\n",
"### State variables\n",
"\n",
"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 \n",
"\n",
"* **velocity**, $v(t)$\n",
"* **altitude**, $h(t)$\n",
"\n",
"Our rocket is burning fuel as it goes, getting lighter and lighter in the process. Therefore we need a third state variable\n",
"\n",
"* **Mass of rocket and remaining fuel**, $m(t)$\n",
"\n",
"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$.\n",
"\n",
"### Control\n",
"\n",
"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}$.\n",
"\n",
"### Dynamics\n",
"\n",
"We need to consider what forces act on our rocket other than thrust. The most obvious perhaps is **gravity**\n",
"\n",
"$$g(h) = g_0 \\left( \\frac{h(0)}{h} \\right)^2$$\n",
"\n",
"which is a function of height. The other key force is **drag**\n",
"\n",
"$$D(h,v) = D_c v^2 exp\\left( -h_c \\left( \\frac{h-h(0)}{h(0)} \\right) \\right)$$\n",
"\n",
"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.\n",
"\n",
"Given all the above, we can link together everything with three equations to 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 $c$ is a constant measuring the impulse of the rocket fuel."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optimization problem\n",
"\n",
"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.\n",
"\n",
"### Dealing with time\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. We can now write out objective as\n",
"$$\n",
"\\max h(t_f)\n",
"$$\n",
"\n",
"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.\n",
"$$\n",
"f(t+\\Delta t) - f(t) = \\Delta t v(t)\n",
"$$\n",
"\n",
"We will use the more accurate but more complicated [trapezoidal rule](http://en.wikipedia.org/wiki/Trapezoidal_rule) (other options are possible!). In particular, in the context our problem we have that\n",
"$$\n",
"f(t+\\Delta t) - f(t) = \\frac{ v(t + \\Delta t) + v(t) }{ 0.5 \\Delta t }\n",
"$$\n",
"\n",
"The three equations that relate the states, controls, and forces correspond to **constraints** in the optimization problem.\n",
"\n",
"### Data\n",
"\n",
"For the purposes of this example we will use dimension free (scaled) parameters. In particular:\n",
"\n",
"* $T_{max} = 3.5 g_0 m_0$\n",
"* $D_c = 0.5 v_c \\frac{m_0}{g_0}$\n",
"* $c = 0.5 (g_0 h_0 )^{0.5}$\n",
"* So $h_0 = m_0 = g_0 = 1$\n",
"* and we choose $h_c = 500$, $m_c = 0.6$, $v_c = 620$\n",
"\n",
"Finally we will use $n$ time steps.\n",
"\n",
"### Initial Solution\n",
"\n",
"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:\n",
"\n",
"* $t_f = n \\Delta t = 1$\n",
"* $h(t) = 1$\n",
"* $v(t) = \\frac{t}{t_f}(1 - \\frac{t}{t_f})$\n",
"* $m(t) = (m_f - m_0)\\frac{t}{t_f} + m_0$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Building the model"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Load JuMP, and a nonlinear solver Ipopt\n",
"# Other options we have for this problem class\n",
"# is KNITRO.jl, a commercial solver.\n",
"using JuMP, Ipopt"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/latex": [
"$$ \\begin{alignat*}{1}\\min\\quad & 0\\\\\n",
"\\text{Subject to} \\quad\\end{alignat*}\n",
" $$"
],
"text/plain": [
"Feasibility problem with:\n",
" * 0 linear constraints\n",
" * 0 variables\n",
"Solver set to Ipopt"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Create JuMP model, using Ipopt as the solver\n",
"#mod = Model(solver=IpoptSolver(print_level=0))\n",
"mod = Model(solver=IpoptSolver())"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"800"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 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",
"# Time steps\n",
"n = 800"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Set up variables and constraints\n",
"\n",
"# Time step with initial guess\n",
"@defVar(mod, Δt ≥ 0, start = 1/n)\n",
"# Store a useful subexpression, \"time of flight\"\n",
"@defNLExpr(t_f, Δt*n)\n",
"\n",
"# State variables\n",
"@defVar(mod, v[0:n] ≥ 0) # Velocity\n",
"@defVar(mod, h[0:n] ≥ h_0) # Height\n",
"@defVar(mod, m_f ≤ m[0:n] ≤ m_0) # Mass\n",
"\n",
"# Control: thrust\n",
"@defVar(mod, 0 ≤ T[0:n] ≤ T_max)\n",
"\n",
"# Provide starting solution\n",
"# Could have done this at same time as @defVar\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",
"# Objective: maximize altitude at end of time of flight\n",
"@setObjective(mod, Max, h[n])\n",
"\n",
"# Initial conditions\n",
"@addConstraint(mod, v[0] == v_0)\n",
"@addConstraint(mod, h[0] == h_0)\n",
"@addConstraint(mod, m[0] == m_0)\n",
"@addConstraint(mod, m[n] == m_f)\n",
"\n",
"# Forces\n",
"# We'll define these as expressions too\n",
"# Wherever they appear, they will effectively be\n",
"# replaced by these longer expressions, keeps them\n",
"# nice and clean\n",
"\n",
"# Drag(h,v) = Dc v^2 exp( -hc * (h - h0) / h0 )\n",
"@defNLExpr(drag[j=0:n], D_c*(v[j]^2)*exp(-h_c*(h[j]-h_0)/h_0))\n",
"# Grav(h) = g0 * (h0 / h)^2\n",
"@defNLExpr(grav[j=0:n], g_0*(h_0/h[j])^2)\n",
"\n",
"# Dynamics\n",
"for j in 1:n\n",
" # h' = v\n",
" # Rectangular integration\n",
" # @addNLConstraint(mod, h[j] == h[j-1] + Δt*v[j-1])\n",
" # Trapezoidal integration\n",
" @addNLConstraint(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",
" # @addNLConstraint(mod, v[j] == v[j-1] + Δt*(\n",
" # (T[j-1] - drag[j-1])/m[j-1] - grav[j-1]))\n",
" # Trapezoidal integration\n",
" @addNLConstraint(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",
" # @addNLConstraint(mod, m[j] == m[j-1] - Δt*T[j-1]/c)\n",
" # Trapezoidal integration\n",
" @addNLConstraint(mod,\n",
" m[j] == m[j-1] - 0.5*Δt*(T[j] + T[j-1])/c)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving...\n",
"This is Ipopt version 3.12.1, running with linear solver mumps.\n",
"NOTE: Other linear solvers might be more efficient (see Ipopt documentation).\n",
"\n",
"Number of nonzeros in equality constraint Jacobian...: 15204\n",
"Number of nonzeros in inequality constraint Jacobian.: 0\n",
"Number of nonzeros in Lagrangian Hessian.............: 48800\n",
"\n",
"Total number of variables............................: 3205\n",
" variables with only lower bounds: 1603\n",
" variables with lower and upper bounds: 1602\n",
" variables with only upper bounds: 0\n",
"Total number of equality constraints.................: 2404\n",
"Total number of inequality constraints...............: 0\n",
" inequality constraints with only lower bounds: 0\n",
" inequality constraints with lower and upper bounds: 0\n",
" inequality constraints with only upper bounds: 0\n",
"\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 0 -1.0100000e+00 3.50e-02 1.83e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n",
" 1 -1.0476357e+00 3.10e-03 7.60e+03 -1.0 2.38e-01 2.0 9.18e-01 9.67e-01f 1\n",
" 2 -1.0330218e+00 1.03e-03 8.61e+04 -1.0 2.31e+00 - 3.71e-01 7.53e-01h 1\n",
" 3 -1.0108102e+00 3.89e-04 1.79e+05 -1.0 1.27e+00 - 1.00e+00 6.40e-01h 1\n",
" 4 -1.0058561e+00 1.58e-04 2.97e+05 -1.0 6.84e-01 - 1.00e+00 9.69e-01h 1\n",
" 5 -1.0078548e+00 9.05e-05 2.16e+05 -1.0 8.08e-01 - 1.00e+00 6.33e-01f 1\n",
" 6 -1.0079653e+00 1.75e-05 4.11e+04 -1.0 4.89e-01 - 1.00e+00 1.00e+00f 1\n",
" 7 -1.0078191e+00 1.38e-06 1.02e+02 -1.0 1.70e-01 - 1.00e+00 1.00e+00f 1\n",
" 8 -1.0078115e+00 3.32e-10 2.15e+02 -1.7 3.11e-03 - 1.00e+00 1.00e+00h 1\n",
" 9 -1.0078157e+00 2.53e-10 4.46e-01 -1.7 1.20e-02 - 1.00e+00 1.00e+00h 1\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 10 -1.0078162e+00 6.07e-12 1.15e+01 -3.8 4.17e-04 - 1.00e+00 1.00e+00h 1\n",
" 11 -1.0078877e+00 1.06e-07 5.24e-03 -3.8 5.52e-02 - 1.00e+00 1.00e+00f 1\n",
" 12 -1.0078881e+00 2.27e-13 1.65e-06 -3.8 3.09e-04 - 1.00e+00 1.00e+00h 1\n",
" 13 -1.0079601e+00 1.06e-07 2.95e+02 -5.7 5.49e-02 - 9.80e-01 1.00e+00f 1\n",
" 14 -1.0095526e+00 8.48e-05 1.22e+00 -5.7 2.23e+00 - 1.00e+00 9.65e-01h 1\n",
" 15 -1.0112897e+00 5.88e-05 2.95e+00 -5.7 1.78e+00 - 1.00e+00 9.90e-01f 1\n",
" 16 -1.0111071e+00 2.22e-06 2.07e-02 -5.7 7.35e-01 - 1.00e+00 1.00e+00f 1\n",
" 17 -1.0111190e+00 6.83e-09 6.56e-05 -5.7 5.17e-02 - 1.00e+00 1.00e+00h 1\n",
" 18 -1.0111191e+00 4.12e-14 1.20e-09 -5.7 3.56e-04 - 1.00e+00 1.00e+00h 1\n",
" 19 -1.0122683e+00 2.30e-05 2.99e+01 -8.6 8.98e-01 - 7.15e-01 8.77e-01f 1\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 20 -1.0127029e+00 1.51e-05 8.14e+00 -8.6 8.55e-01 - 7.50e-01 7.80e-01h 1\n",
" 21 -1.0128032e+00 1.01e-05 3.13e+00 -8.6 1.26e+00 - 6.59e-01 7.31e-01h 1\n",
" 22 -1.0128269e+00 5.32e-06 1.22e+00 -8.6 1.40e+00 - 6.46e-01 7.33e-01h 1\n",
" 23 -1.0128326e+00 2.58e-06 3.96e-01 -8.6 1.38e+00 - 6.98e-01 7.71e-01f 1\n",
" 24 -1.0128339e+00 1.13e-06 1.87e-02 -8.6 1.24e+00 - 9.40e-01 8.79e-01f 1\n",
" 25 -1.0128341e+00 2.77e-07 5.12e-05 -8.6 9.73e-01 - 1.00e+00 1.00e+00f 1\n",
" 26 -1.0128341e+00 3.97e-09 1.13e-06 -8.6 2.34e-01 - 1.00e+00 1.00e+00h 1\n",
" 27 -1.0128341e+00 1.03e-10 7.14e-09 -8.6 3.13e-02 - 1.00e+00 1.00e+00h 1\n",
"\n",
"Number of Iterations....: 27\n",
"\n",
" (scaled) (unscaled)\n",
"Objective...............: -1.0128340615485458e+00 -1.0128340615485458e+00\n",
"Dual infeasibility......: 7.1350593378980852e-09 7.1350593378980852e-09\n",
"Constraint violation....: 1.0326103860869296e-10 1.0326103860869296e-10\n",
"Complementarity.........: 2.5081820874400327e-09 2.5081820874400327e-09\n",
"Overall NLP error.......: 7.1350593378980852e-09 7.1350593378980852e-09\n",
"\n",
"\n",
"Number of objective function evaluations = 28\n",
"Number of objective gradient evaluations = 28\n",
"Number of equality constraint evaluations = 28\n",
"Number of inequality constraint evaluations = 0\n",
"Number of equality constraint Jacobian evaluations = 28\n",
"Number of inequality constraint Jacobian evaluations = 0\n",
"Number of Lagrangian Hessian evaluations = 27\n",
"Total CPU secs in IPOPT (w/o function evaluations) = 0.259\n",
"Total CPU secs in NLP function evaluations = 0.265\n",
"\n",
"EXIT: Optimal Solution Found.\n",
"Solver status: Optimal\n",
"Max height: 1.0128340615485458\n"
]
}
],
"source": [
"# Solve for the control and state\n",
"println(\"Solving...\")\n",
"status = solve(mod)\n",
"\n",
"# Display results\n",
"println(\"Solver status: \", status)\n",
"println(\"Max height: \", getObjectiveValue(mod))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Can visualize the state and control variables\n",
"using Gadfly"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"SVG(152.39999999999998,152.39999999999998,IOBuffer(data=Uint8[...], readable=true, writable=true, seekable=true, append=false, size=54579, maxsize=Inf, ptr=54580, mark=-1),nothing,\"fig-2a34be1329b044df9921e3cf270e22e3\",0,[],Dict{Type{T<:Top},Union(Nothing,Property{P<:PropertyPrimitive})}(),[ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(21.97833333333333,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(21.97833333333333,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-42\",ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(20.305,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(71.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(20.305,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-18\",ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(94.83166666666665,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(5.0,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(94.83166666666665,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(56.914999999999985,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-30\",ClipPrimitive{Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}}([Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(92.48499999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(81.19999999999999,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(147.39999999999998,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0)),Point{Measure{MeasureNil,MeasureNil},Measure{MeasureNil,MeasureNil}}(Measure{MeasureNil,MeasureNil}(92.48499999999999,MeasureNil(),MeasureNil(),0.0,0.0),Measure{MeasureNil,MeasureNil}(133.11499999999998,MeasureNil(),MeasureNil(),0.0,0.0))])=>\"fig-2a34be1329b044df9921e3cf270e22e3-element-6\"],Set{String}({}),true,false,nothing,true,\"fig-2a34be1329b044df9921e3cf270e22e3-element-48\",false,48,String[\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\",\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\",\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\",\"/Users/idunning/.julia/v0.3/Gadfly/src/gadfly.js\"],(String,String)[(\"Snap.svg\",\"Snap\"),(\"Gadfly\",\"Gadfly\"),(\"Gadfly\",\"Gadfly\"),(\"Gadfly\",\"Gadfly\"),(\"Gadfly\",\"Gadfly\")],String[\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-5\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\",\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-17\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\",\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-29\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\",\"fig.select(\\\"#fig-2a34be1329b044df9921e3cf270e22e3-element-41\\\")\\n .mouseenter(Gadfly.plot_mouseover)\\n.mouseleave(Gadfly.plot_mouseout)\\n.mousewheel(Gadfly.guide_background_scroll)\\n.drag(Gadfly.guide_background_drag_onmove,\\n Gadfly.guide_background_drag_onstart,\\n Gadfly.guide_background_drag_onend)\\n;\"],false,:none)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: groupby(xs, keyfunc) should be groupby(keyfunc, xs)\n"
]
}
],
"source": [
"h_plot = plot(x=(0:n)*getValue(Δt),y=getValue(h)[1:n], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Altitude\"))\n",
"m_plot = plot(x=(0:n)*getValue(Δt),y=getValue(m)[1:n], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Mass\"))\n",
"v_plot = plot(x=(0:n)*getValue(Δt),y=getValue(v)[1:n], Geom.line,\n",
" Guide.xlabel(\"Time (s)\"), Guide.ylabel(\"Velocity\"))\n",
"T_plot = plot(x=(0:n)*getValue(Δt),y=getValue(T)[1:n], 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)))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.3.7",
"language": "julia",
"name": "julia 0.3"
},
"language_info": {
"name": "julia",
"version": "0.3.7"
}
},
"nbformat": 4,
"nbformat_minor": 0
}