Example about climbing (increase altitude)
Reference: [5]
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;
sos1 = interp1(hTab,sosTab,h*0.3048,'pchip')./0.3048;
Mach = v/sos1;
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';
% Strating guess of fpa is in confilct with the boundary conditions, but
% that's ok. (It will give a warning, which we suppress.)
warns = warning('off', 'tomSym:x0OutOfBounds');
ezsolve(t_f,{cbox,bnd,equ},guess,options);
% Restore warning
warning(warns);
Problem type appears to be: lpcon
Auto-scaling
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - Tomlab Optimization Inc. Development license 999001. Valid to 2011-02-05
=====================================================================================
Problem: --- 1: Minimum Time to Climb (English) f_k 318.348640856079270000
sum(|constr|) 0.000000000036130174
f(x_k) + sum(|constr|) 318.348640856115420000
f(x_0) 300.000000000000000000
Solver: snopt. EXIT=0. INFORM=3.
SNOPT 7.2-5 NLP code
Requested accuracy could not be achieved
FuncEv 1 ConstrEv 96 ConJacEv 95 Iter 51 MinorIter 2042
CPU time: 4.703125 sec. Elapsed time: 2.969000 sec.
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');