PROPT Minimum Climb Time (English Units): Difference between revisions

From TomWiki
Jump to navigationJump to search
(Created page with "{{Part Of Manual|title=the PROPT Manual|link=PROPT Manual}} ==Problem description== Example about climbing (increase altitude) <source lang="matlab"> % Copyright ...")
 
No edit summary
Line 1: Line 1:
{{Part Of Manual|title=the PROPT Manual|link=[[PROPT |PROPT Manual]]}}
{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}


==Problem description==
==Problem description==
Line 174: Line 174:
title('Minimum Time to Climb (English Units) control');
title('Minimum Time to Climb (English Units) control');
</source>
</source>
[[File:minimumClimbEng_01.png]]

Revision as of 14:24, 2 November 2011

Notice.png

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

Problem description

Example about climbing (increase altitude)

% Copyright (c) 2007-2008 by Tomlab Optimization Inc.

Problem setup

alt0     = 0;
altf     = 65600;
speed0   = 424.26;
speedf   = 968.148;
fpa0     = 0;
fpaf     = 0;
mass0    = 42000/32.208;

altmin   = 0;
altmax   = 69000;
speedmin = 10;
speedmax = 3000;
fpamin   = -40*pi/180;
fpamax   = -fpamin;
massmin  = 50/32.208;

toms t t_f
p = tomPhase('p',t,0,t_f,50);
setPhase(p);

% Altitude, speed, flight path angle, mass
tomStates h v fpa m

% Angle of Attack
tomControls aalpha

guess = {
    t_f == 300;
    icollocate({
    h == alt0 + t/t_f*(altf-alt0);
    v == speed0 + t/t_f*(speedf-speed0);
    fpa == 10*pi/180;
    m == mass0;
    })};

cbox = {
    100 <= t_f <= 800;
    collocate(-pi*20/180 <= aalpha <= pi/20*180)
    icollocate(altmin <= h <= altmax)
    icollocate(speedmin <= v <= speedmax)
    icollocate(fpamin <= fpa <= fpamax)
    icollocate(massmin <= m <= mass0)
    };

bnd = {
    initial(h) == alt0;
    initial(v) == speed0;
    initial(fpa) == fpa0;
    initial(m) == mass0;
    final(h) == altf;
    final(v) == speedf;
    final(fpa) == fpaf;
    };

% US1976 data
hTab = (-2000:2000:86000);
rhoTab = [1.478 1.225 1.007 0.8193 0.6601 0.5258 0.4135 0.3119 ...
    0.2279 0.1665 0.1216 0.08891 0.06451 0.04694 0.03426 0.02508 ...
    0.01841 0.01355 0.009887 0.007257 0.005366 0.003995 0.002995 ...
    0.002259 0.001714 0.001317 0.001027 0.0008055 0.0006389 0.0005044 ...
    0.0003962 0.0003096 0.0002407 0.000186 0.0001429 0.0001091 ...
    8.281e-005 6.236e-005 4.637e-005 3.43e-005 2.523e-005 1.845e-005 ...
    1.341e-005 9.69e-006 6.955e-006];
sosTab = [347.9 340.3 332.5 324.6 316.5 308.1 299.5 295.1 295.1 ...
    295.1 295.1 295.1 296.4 297.7 299.1 300.4 301.7 303 306.5 310.1 ...
    313.7 317.2 320.7 324.1 327.5 329.8 329.8 328.8 325.4 322 318.6 ...
    315.1 311.5 308 304.4 300.7 297.1 293.4 290.7 288 285.3 282.5 ...
    279.7 276.9 274.1];

Mtab   = [0; 0.2; 0.4; 0.6; 0.8; 1; 1.2; 1.4; 1.6; 1.8];
alttab = [0 5000 10000 15000 20000 25000 30000 40000 50000 70000];
Ttab = 1000*[24.2 24.0  20.3 17.3 14.5 12.2 10.2 5.7 3.4 0.1;
    28.0 24.6 21.1 18.1 15.2 12.8 10.7 6.5 3.9 0.2;
    28.3 25.2 21.9 18.7 15.9 13.4 11.2 7.3 4.4 0.4;
    30.8 27.2 23.8 20.5 17.3 14.7 12.3 8.1 4.9 0.8;
    34.5 30.3 26.6 23.2 19.8 16.8 14.1 9.4 5.6 1.1;
    37.9 34.3 30.4 26.8 23.3 19.8 16.8 11.2 6.8 1.4;
    36.1 38.0 34.9 31.3 27.3 23.6 20.1 13.4 8.3 1.7;
    36.1 36.6 38.5 36.1 31.6 28.1 24.2 16.2 10.0 2.2;
    36.1 35.2 42.1 38.7 35.7 32.0 28.1 19.3 11.9 2.9;
    36.1 33.8 45.7 41.3 39.8 34.6 31.1 21.7 13.3 3.1];
M2     = [0 0.4 0.8 0.9 1.0 1.2 1.4 1.6 1.8];
Clalphatab = [3.44 3.44 3.44 3.58 4.44 3.44 3.01 2.86 2.44];
CD0tab    = [0.013 0.013 0.013 0.014 0.031 0.041 0.039 0.036 0.035];
etatab    = [0.54 0.54 0.54 0.75 0.79 0.78 0.89 0.93 0.93];
M         = Mtab;
alt       = alttab;
Re        = 20902900;
mmu        = 0.14076539e17;
S         = 530;
g0        = 32.208;
ISP       = 1600;
H         = 23800;
rho0      = 0.002378;

rho = interp1(hTab,rhoTab,h*0.3048,'pchip')*0.001941;
sosI = interp1(hTab,sosTab,h*0.3048,'pchip')./0.3048;
Mach = v/sosI;

CD0       = interp1(M2,CD0tab,Mach,'pchip');
Clalpha   = interp1(M2,Clalphatab,Mach,'pchip');
eta       = interp1(M2,etatab,Mach,'pchip');

T = interp2(alttab, Mtab, Ttab, h, Mach, 'spline');
CD      = CD0 + eta.*Clalpha.*aalpha.^2;
CL      = Clalpha.*aalpha;
dynpres = 0.5.*rho.*v.^2;
D       = dynpres.*S.*CD;
L       = dynpres.*S.*CL;

equ = collocate({
    dot(h) == v.*sin(fpa);
    dot(v) == ((T.*cos(aalpha)-D)./m - mmu.*sin(fpa)./(Re+h).^2);
    dot(fpa) == (T.*sin(aalpha)+L)./(m.*v)+cos(fpa).*(v./(Re+h)-mmu./(v.*(Re+h).^2));
    dot(m) == -T./(g0.*ISP);
    });

options = struct;
options.name = 'Minimum Time to Climb (English)';
options.scale = 'auto';

Solve the problem

ezsolve(t_f,{cbox,bnd,equ},guess,options);
Problem type appears to be: lpcon
Auto-scaling
Time for symbolic processing: 0.89551 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Minimum Time to Climb (English)  f_k     318.348640855936940000
                                         sum(|constr|)      0.000000000187118413
                                f(x_k) + sum(|constr|)    318.348640856124060000
                                                f(x_0)    300.000000000000000000

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

FuncEv    1 ConstrEv   91 ConJacEv   90 Iter   49 MinorIter 2246
CPU time: 2.168414 sec. Elapsed time: 1.359000 sec. 

Plot result

subplot(2,1,1)
ezplot([1e-3*h, 1e-1*v, fpa*180/pi, 1e-1*m])
legend('1e-3*h', '1e-1*v', 'fpa*180/pi', '1e-1*m');
title('Minimum Time to Climb (English Units) state variables');

subplot(2,1,2)
ezplot(aalpha);
title('Minimum Time to Climb (English Units) control');

MinimumClimbEng 01.png