PROPT Goddard Rocket, Maximum Ascent, Final time free, Singular solution

This page is part of the PROPT Manual. See PROPT Manual.

Example 7.2 (i) from the paper: H. Maurer, "Numerical solution of singular control problems using multiple shooting techniques", Journal of Optimization Theory and Applications, Vol. 18, No. 2, 1976, pp. 235-257

L.G. van Willigenburg, W.L. de Koning

Copyright (c) 2007-2009 by Tomlab Optimization Inc.

Problem setup

toms t t_f

% Parameters
aalpha = 0.01227; bbeta = 0.145e-3;

c  = 2060;    g0 = 9.81;
r0 = 6.371e6; r02=r0*r0;
m0 = 214.839; mf = 67.9833;
Fm = 9.525515;

Solve the problem, using a successively larger number collocation points

for n=[20 40 60]

    p = tomPhase('p', t, 0, t_f, n);
    setPhase(p);
    tomStates h v m
    tomControls F

    % Initial guess
    if n==20
        x0 = {t_f==250
            icollocate({v == 0; h == 0
            m == m0})
            collocate(F == Fm)};
    else
        x0 = {t_f == tfopt
            icollocate({v == vopt; h == hopt
            m == mopt})
            collocate(F == Fopt)};
    end

    % Box constraints
    cbox = {100 <= t_f <= 300
        icollocate({0 <= v; 0 <= h
        mf <= m <= m0
        0 <= F <= Fm})};

    % Boundary constraints
    cbnd = {initial({v == 0; h == 0; m == m0})
        final({v==0; m == mf})};

    D = aalpha*v.^2.*exp(-bbeta*h);
    g = g0; % or g0*r02./(r0+h).^2;

    % ODEs and path constraints
    ceq = collocate({dot(h) == v
        m*dot(v) == F*c-D-g*m
        dot(m) == -F});

    % Objective
    objective = -1e-4*final(h);

Solve the problem

    options = struct;
    options.name = 'Goddard Rocket 1';
    options.Prob.SOL.optPar(30) = 30000;
    solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

    % Optimal v and more to use as starting guess
    vopt = subs(v, solution);
    hopt = subs(h, solution);
    mopt = subs(m, solution);
    Fopt = subs(F, solution);
    tfopt = subs(t_f, solution);

Problem type appears to be: lpcon
Time for symbolic processing: 0.16423 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Goddard Rocket 1               f_k     -15.580049356479096000
                                       sum(|constr|)      0.000028629544147830
                              f(x_k) + sum(|constr|)    -15.580020726934949000
                                              f(x_0)      0.000000000000000000

Solver: snopt.  EXIT=0.  INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv    1 ConstrEv  217 ConJacEv  216 Iter   45 MinorIter 1272
CPU time: 0.140401 sec. Elapsed time: 0.151000 sec. 

Problem type appears to be: lpcon
Time for symbolic processing: 0.15968 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Goddard Rocket 1               f_k     -15.718139470225974000
                                       sum(|constr|)      0.043626439767554738
                              f(x_k) + sum(|constr|)    -15.674513030458419000
                                              f(x_0)    -15.580049356479037000

Solver: snopt.  EXIT=0.  INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv    1 ConstrEv  235 ConJacEv  234 Iter   34 MinorIter  915
CPU time: 0.280802 sec. Elapsed time: 0.290000 sec. 

Problem type appears to be: lpcon
Time for symbolic processing: 0.16035 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Goddard Rocket 1               f_k     -15.731752507215452000
                                       sum(|constr|)      0.000170868961459850
                              f(x_k) + sum(|constr|)    -15.731581638253992000
                                              f(x_0)    -15.718139470225950000

Solver: snopt.  EXIT=0.  INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv    1 ConstrEv  234 ConJacEv  233 Iter   53 MinorIter 3497
CPU time: 0.904806 sec. Elapsed time: 0.917000 sec. 

end

t = subs(collocate(t),solution);
v = subs(collocate(vopt),solution);
h = subs(collocate(hopt),solution);
m = subs(collocate(mopt),solution);
F = subs(collocate(Fopt),solution);

Plot result

subplot(2,1,1)
plot(t,v/1e3,'*-',t,h/1e5,'*-',t,m/1e2,'*-');
legend('v','h','m');
title('Goddard Rocket state variables');

subplot(2,1,2)
plot(t,F,'+-');
legend('F');
title('Goddard Rocket control');