Goddard rocket

goddard.jl

#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# goddard.jl
#
# Goddard Rocket control problem from:
#  Benchmarking Optimization Software with COPS 3.0
#  http://www.mcs.anl.gov/~more/cops/cops3.pdf
#############################################################################

using JuMP, Ipopt

# Constants
h_0 = 1    # Initial height
v_0 = 0    # Initial velocity
m_0 = 1    # Initial mass
g_0 = 1    # Gravity at the surface

# Parameters
# Not very interpretable, due to everything being dimensionless
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

n   = 400   # Time steps
Δt  = 1/n   # Time step
t_f = Δt*n  # Time of flight

# Create JuMP model, using Ipopt as the solver
mod = Model(solver=IpoptSolver())

# 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)

# 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
# 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)   = go * (h0 / h)^2
@defNLExpr(grav[j=0:n], g_0*(h_0/h[j])^2)

# Dynamics - using trapezoidal rule
for j in 1:n
    # h' = v
    @addNLConstraint(mod, h[j] == h[j-1] + 0.5*Δt*(v[j] + v[j-1]))

    # v' = (T-D(h,v))/m - g(h)
    @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
    @addNLConstraint(mod, m[j] == m[j-1] - 0.5*Δt*(T[j] + T[j-1])/c)
end

# Provide starting solution
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

# Solve for the control and state
println("Solving...")
status = solve(mod)

# Display results
println("Solver status: ", status)
println("Max height: ", getObjectiveValue(mod))

rocket.dat

param v_0 := 0;

# Normalization of the equations allows g_0 = h_0 = m_0 = 1

param g_0 := 1.0;
param h_0 := 1.0;
param m_0 := 1.0;

param T_c := 3.5;
param h_c := 500;
param v_c := 620;
param m_c := 0.6;

# Starting values

let {k in 0..nh} h[k] := 1;;
let {k in 0..nh} v[k] := (k/nh)*(1 - (k/nh));
let {k in 0..nh} m[k] := (m_f - m_0)*(k/nh) + m_0;
let {k in 0..nh} T[k] := T_max/2;
let step := 1/nh;

rocket.mod

# Goddard Rocket Problem
# Trapezoidal formulation
# COPS 2.0 - September 2000
# COPS 3.0 - November 2002
# COPS 3.1 - March 2004

param h_0;     # Initial height
param v_0;     # Initial velocity
param m_0;     # Initial mass
param g_0;     # Gravity at the surface

# Parameters for the model.

param T_c;
param h_c;
param v_c;
param m_c;

# Derived parameters.

param c := 0.5*sqrt(g_0*h_0);
param m_f := m_c*m_0;
param D_c := 0.5*v_c*(m_0/g_0);
param T_max := T_c*(m_0*g_0);

param nh;      # Number of intervals in mesh

# Height, velocity, mass and thrust of rocket.

var h {0..nh};
var v {0..nh};
var m {0..nh};
var T {0..nh};
var step;
var t_f = step*nh;

# Drag function.

var D {i in 0..nh} = D_c*(v[i]^2)*exp(-h_c*(h[i]-h_0)/h_0);

# Gravity function.

var g{i in 0..nh} = g_0*(h_0/h[i])^2;

maximize final_velocity: h[nh];

subject to step_eqn: step >= 0;

subject to v_bounds {j in 0..nh}: v[j] >= 0.0;
subject to h_bounds {j in 0..nh}: h[j] >= h_0;
subject to m_bounds {j in 0..nh}: m_f <= m[j] <= m_0;
subject to T_bounds {j in 0..nh}: 0.0 <= T[j] <= T_max;

subject to h_eqn {j in 1..nh}:
   h[j] = h[j-1] + .5*step*(v[j] + v[j-1]);

subject to v_eqn {j in 1..nh}:
   v[j] = v[j-1] + .5*step*((T[j] - D[j] - m[j]*g[j])/m[j]
                 +          (T[j-1] - D[j-1] - m[j-1]*g[j-1])/m[j-1]);

subject to m_eqn {j in 1..nh}:
   m[j] = m[j-1] - .5*step*(T[j] + T[j-1])/c;

# Boundary Conditions

subject to h_ic : h[0] = h_0;
subject to v_ic : v[0] = v_0;
subject to m_ic : m[0] = m_0;
subject to m_fc : m[nh] = m_f;