Rocket : Dynamic Optimization of a Rocket

Reference

• Neculai Andrei, Nonlinear Optimization Applications Using the GAMS Technology,Springer Optimization and Its Applications, Model Rocket (5.30) in chapter Applications of Mechanical Engineering , 2013

Category : GAMS NOA library

Mainfile : rocket.gms

``````\$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
\$onecho >bench.opt
solvers conopt knitro minos
\$offecho

rocket.optfile=1;
rocket.iterlim=80000;

option nlp=bench;
\$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
*------------------------------------------------------------------
``````