$onText
    Goddard Rocket.
    Maximize the final altitude of a vertically launched rocket, using the
    thrust as a control and given the initial mass, the fuel mass, and the
    drag characteristics of the rocket.

    This model is from the COPS benchmarking suite.
    See http://www-unix.mcs.anl.gov/~more/cops/.

    The number of discretization points can be specified using the
    command line parameter.
    COPS performance tests have been reported for nh = 50,
    100, 200, 400

    References:
  * Dolan, E D, and More, J J, Benchmarking Optimization
    Software with COPS. Tech. rep., Mathematics and Computer
    Science Division, 2000.
  * Bryson, A E, Dynamic Optimization. Addison Wesley, 1999.
$offText

$if     set n  $set nh %n%
$if not set nh $set nh 1000

set h intervals / h0 * h%nh%/

scalars
   h_0  Initial height          / 1 /
   v_0  Initial velocity        / 0 /
   m_0  Initial mass            / 1 /
   g_0  Gravity at the surface  / 1 /
   nh   Number of intervals in mesh / %nh% /
   r_c  Thrust constant         /3.5/
   v_c  / 620 /
   h_c  / 500 /
   m_c  / 0.6 /
   D_c
   m_f  final mass
   c ;

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

variable   final_velocity

positive
variables  step   step size
           v(h)   velocity
           ht(h)  height
           g(h)   gravity
           m(h)   mass
           r(h)   thrust
           d(h)   drag     ;

* Bounds:
ht.lo(h) = h_0;

r.lo(h)  = 0.0;
r.up(h)  = r_c*(m_0*g_0);

m.lo(h)  = m_f;
m.up(h)  = m_0;

* Initial values:
ht.l(h) = 1;
v.l(h)  = ((ord(h)-1)/nh)*(1 - ((ord(h)-1)/nh));
m.l(h)  = (m_f - m_0)*((ord(h)-1)/nh) + m_0;
r.l(h)  = r.up(h)/2;
step.l  = 1/nh;

d.l(h) = D_c*sqr(v.l(h))*exp(-h_c*(ht.l(h)-h_0)/h_0);
g.l(h) = g_0*sqr(h_0/ht.l(h));

* Fixed values for variables:
ht.fx('h0')   = h_0;
v.fx('h0')    = v_0;
m.fx('h0')    = m_0;
m.fx('h%nh%') = m_f;

equations  df(h)  Drag function
           gf(h)  Gravity function
           obj
           h_eqn(h), v_eqn(h), m_eqn(h);

obj.. final_velocity =e= ht('h%nh%');

df(h).. d(h) =e= D_c*sqr(v(h))*exp(-h_c*(ht(h)-h_0)/h_0);

gf(h).. g(h) =e= g_0*sqr(h_0/ht(h));

h_eqn(h-1).. ht(h) =e= ht(h-1) + .5*step*(v(h) + v(h-1));

m_eqn(h-1).. m(h) =e= m(h-1) - .5*step*(r(h) + r(h-1))/c;

v_eqn(h-1).. v(h) =e= v(h-1)
                    + .5*step*((r(h)   - D(h)   - m(h)  *g(h))  /m(h)
                              +(r(h-1) - D(h-1) - m(h-1)*g(h-1))/m(h-1));

model rocket /all/;

$ifThenI x%mode%==xbook
rocket.iterlim=80000;
$endIf

solve rocket using nlp maximizing final_velocity;

$ifThenI x%mode%==xbook
 file res1 /m9.dat/;
 put res1

 loop(h, put r.l(h):10:7, put/)
$endIf

* End rocket

*------------------ Numerical Experiments I have obtained ----------------
*                             January 15, 2011
*
*
*                                Variant 1:
* nh=500
*       m=2502,  n=3007
* CONOPT3:
* 1997 iterations,  22.713 seconds
* vfo=1.012836666
* KNITRO:
* 70 major iterations,  91 minor iterations,  171 functions evaluations
* 5.89 seconds
* vfo=1.01262977
* MINOS:
* iter=43, fev=2781, 0.230 seconds
* vfo=1.102603
*
*                                Variant 2:
* nh=1000
*        m=5003,   n=6008
* CONOPT3:
* 2700 iterations,  49.801 seconds
* vfo=1.012835932
* KNITRO:
* 71 major iterations, 93 minor iterations, 164 functions evaluations
* 19.327 seconds
* vfo=1.0123973821
* MINOS:
* 35 iterations,  38.846 seconds
* vfo=1.012243
*------------------------------------------------------------------