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');