-

Advanced_control_of_aircraft_spacecraft_Tewari/2.5.m

@@ -0,0 +1,19 @@
+% Calling program for solving the 2PBVP of a second-order
+% system d^2y/dt^2=f(t,y,dy/dt),
+% by collocation method of MATLAB's intrinsic code 'bvp4c.m'.
+% Requires 'ode2ord.m' and 'bc2ord.m'.
+% (c) 2009 Ashish Tewari
+% Collocation solution with 5 points with t0=0, tf=1:
+solinit = bvpinit(linspace(0,1,5),[1 0]); % Initial guess
+%
+(y=1, dy/dt=0)
+sol = bvp4c(@ode2ord,@bc2ord,solinit);% Solution by bvp4c.m
+t = linspace(0,1); % Time vector
+x = deval(sol,t); % Solution state vector, [y, dy/dt]'
+% Program for specifying boundary conditions for the 2PBVP.
+% (To be called by 'bvp4c.m')
+function res=bc2ord(ya,yb)
+% Boundary conditions follow:
+res=[ya(1)-1
+% y(t0)=1
+yb(1)-2]; % y(tf)=2

Advanced_control_of_aircraft_spacecraft_Tewari/3.1.m

@@ -0,0 +1,35 @@
+%Calling program for solving the 2PBVP for optimal aircraft
+%navigation by collocation method using MATLABs
+%intrinsic code bvp4c.m.
+%Requires airnavode.m and airnavbc.m.
+%(c) 2010 Ashish Tewari
+%dy/dx=f(y,x); a<=x<=b
+%y(x=a), y(x=b): Boundary conditions
+%y(1,1)= Latitude (rad.)
+%y(2,1)= Longitude (rad.)
+%y(3,1)= Velocity azimuth (rad.)
+%y(4:6,1)=lambda (Lagrange multipliers vector)
+global del0; del0=51.5*pi/180; %initial latitude
+global lam0;lam0=0; %initial longitude
+global delf; delf=40.4*pi/180; %final latitude
+global lamf;lamf=-73.5*pi/180; %final longitude
+global A0;A0=287.2*pi/180; %initial velocity azimuth
+global vw;vw=20; %wind speed
+global Aw;Aw=45*pi/180; %wind velocity azimuth
+global R; R=0.001; %control cost coefficient
+global vprime; vprime=270; %airspeed (m/s)
+global r; r=6378.14e3+11000; %radius (m)
+global Omega; Omega=7.2921e-5; %Earths rotational rate (rad/s)
+tf=6.12*3600; %terminal time
+
+% Collocation points & initial guess follow:
+solinit = bvpinit(linspace(0,tf,50),[del0 lam0 A0 0 0 0]);
+% 2PBVP solution by collocation method:
+sol = bvp4c(@airnavode,@airnavbc,solinit);
+x = linspace(0,tf);
+y = deval(sol,x);
+% Ground speed (m/s) follows:
+v=sqrt(vprime^2-vw^2*sin(Aw-y(3,:)).^2)+vw*cos(Aw-y(3,:));
+% Normal acceleration input and bank angle profile follow:
+u=-y(6,:)./(2*R*v);
+sigma=atan(u/9.81);

Advanced_control_of_aircraft_spacecraft_Tewari/3.4.m

@@ -0,0 +1,39 @@
+%Calling program for simulating the great circle trajectory for
+%navigation by fourth-order Runge-Kutta solution.
+%Requires airnav.m.
+%(c) 2010 Ashish Tewari
+%y(1,1)= Latitude (rad.)
+%y(2,1)= Longitude (rad.)
+%y(3,1)= Velocity azimuth (rad.)
+global del0; del0=51.5*pi/180; %initial latitude
+global lam0;lam0=0; %initial longitude
+global delf; delf=40.4*pi/180; %final latitude
+global lamf;lamf=-73.5*pi/180; %final longitude
+global vw;vw=20; %wind speed
+global Aw;Aw=45*pi/180; %wind velocity azimuth
+global vprime; vprime=270; %airspeed (m/s)
+global r; r=6378.14e3+11000; %radius (m)
+global Omega; Omega=7.2921e-5; %Earths rotational rate (rad/s)
+tf=6.12*3600; %terminal time
+dtr=pi/180;
+% Computation of the initial velocity azimuth, A_0:
+r0=[cos(del0)*cos(lam0) cos(del0)*sin(lam0) sin(del0)]
+rf=[cos(delf)*cos(lamf) cos(delf)*sin(lamf) sin(delf)]
+n=cross(r0,rf)/norm(cross(r0,rf))
+cosi=dot([0 0 1],n)
+sinA0=cosi/cos(del0);
+A0=asin(sinA0)
+v0=sin(A0)*[-sin(lam0) cos(lam0) 0]+ ...
+cos(A0)*[-sin(del0)*cos(lam0) ...
+-sin(del0)*sin(lam0) cos(del0)];
+h0=cross(r0,v0)/norm(cross(r0,v0))
+if abs(dot(h0,n))<=0.99
+A0=pi-A0
+end
+% 4th order, 5th stage Runge-Kutta integration:
+[x,y]=ode45(@airnav,[0 tf],[del0 lam0 A0]');
+% Ground speed (m/s) follows:
+v=sqrt(vprime^2-vw^2*sin(Aw-y(:,3)).^2)+vw*cos(Aw-y(:,3));
+% Normal acceleration input and bank angle profile follow:
+u=-0.5*Omega^2*r*sin(2*y(:,1)).*sin(y(:,3))-2*Omega*v.*sin(y(:,1));
+sigma=atan(u/9.81);

Advanced_control_of_aircraft_spacecraft_Tewari/3.6.m

@@ -0,0 +1,34 @@
+%Calling program for solving the 2PBVP for optimal aircraft
+%rotational maneuver by collocation method using MATLABs
+%intrinsic code 
+%Requires 'airrotode.m' and 'airrotbc.m'.
+%(c) 2009 Ashish Tewari
+%dy/dx=f(y,x); a<=x<=b
+%y(x=a), y(x=b): Boundary conditions
+%y(1:3,1)=omega (rad/s) (Angular velocity vector)
+%y(4:7,1)=q_1,q_2,q_3,q_4 (Quaternion)
+%y(8:14,1)=lambda (Lagrange multipliers vector)
+global tf; tf=2; % Terminal time (s)
+global J; J=[15000 0 500;0 36000 0;500 0 50000];%(kg-m^2)
+global wi; wi=[1 0.3 -0.02]'; % Initial ang. velocity (rad/s)
+global Qi; Qi=[0 0 0 1]'; % Initial quaternion
+global R; R=inv(J); % Control cost coefficients
+% Collocation points & initial guess follow:
+solinit = bvpinit(linspace(0,tf,5),[zeros(3,1);
+[0.125 0 0 sqrt(1-.125^2)]'; zeros(7,1)]);
+% 2PBVP Solution by collocation method:
+options=bvpset('Nmax',100);
+sol = bvp4c(@airrotode,@airrotbc,solinit,options);
+x = linspace(0,tf); % Time vector (s)
+y = deval(sol,x); % Solution state vector
+plot(x,y(1:3,:)),xlabel('Time (s)'),ylabel('\omega (rad/s)')
+figure
+
+plot(x,y(4:7,:)),xlabel('Time (s)'),ylabel('q,q_4')
+figure
+jx=J(1,1);jy=J(2,2);jz=J(3,3);jxz=J(1,3);
+Jxz=J(1,3);
+D=jx*jz-jxz^2;
+Fu=[jz/D 0 -jxz/D; 0 1/jy 0; -jxz/D 0 jx/D];
+u=-0.5*inv(Rbar)*Fu'*y(8:10,:);
+plot(x,u/1000),xlabel('Time (s)'),ylabel('u (x1000 N-m)')

Advanced_control_of_aircraft_spacecraft_Tewari/4.11.m

@@ -0,0 +1,25 @@
+% Calling program for solving the 2PBVP for
+% flight with modulated fwd and normal
+% acceleration inputs by collocation method
+% using MATLAB's intrinsic code 'bvp4c.m'.
+% Requires 'normballisticode.m' and
+% 'normballisticbc.m'.
+% (c) 2010 Ashish Tewari
+% dy/dx=f(y,x); a<=x<=b
+% y(x=a), y(x=b): Boundary conditions
+% y(1,1)=delta (rad)
+% y(2,1)=r (km)
+% y(3,1)=v (km/s)
+% y(4,1)=phi (rad)
+% y(5:8,1)=Lagrange multipliers
+global R1; R1=1;
+global R2; R2=7;
+global tf; tf=500;
+dtr=pi/180;
+init=[30*dtr 6580 8 5*dtr 0 0 0 0];
+solinit = bvpinit(linspace(0,tf,5),init);
+sol = bvp4c(@normballisticode,@normballisticbc,solinit);
+x = linspace(0,tf);
+y = deval(sol,x);
+u1=-0.5*y(7,:)/R1;
+u2=-0.5*y(8,:)./(R2*y(3,:));

Advanced_control_of_aircraft_spacecraft_Tewari/4.3.m

@@ -0,0 +1,34 @@
+global T;T=[0 1 2 3 4 5 6 7 10 15 20 25 30]';
+global aT;aT=[-0.5 1 3 3.5 3.6 3.8 3.9 4 4.11 4.12 4.13 4.15 4.2
+0.1 -1 -5 -5.5 -5.6 -5.8 -6.0 -5.9 -5.2 -4 -3.5 -3.6 -3.7
+2.5 2 1 -1 -1.5 -2 -1.6 0.4 1.6 2.1 2.7 2.8 3]';
+global kr;kr=.2;
+global kv;kv=.7;
+global M; M = [0; 0.2; 0.3; 0.4; 0.5; 0.6; 0.8; 0.9; 0.95; 1.05; 1.1; 1.2;
+1.6; 2.0; 2.5; 3; 3.8; 5; 10; 99];
+global CD;CD=[0.14264 0.14264 0.142728 0.145008 0.146688 0.152504 ...
+0.169688
+0.185448
+0.20044
+0.309536 0.30512 0.297168 ...
+0.244784 0.207712 0.182912
+0.163816
+0.153904
+0.1482 ...
+0.144952 0.144952];
+global S; S = 0.09; % reference base area, m2
+global tb; tb=20; % burn-time (s)
+global m0; m0=600; % initial mass (kg)
+global mf; mf=80; % burn-out mass (kg)
+[t,x]=ode45('missile_state',[0 20],[1000 2000 3000 100 -50 10 0 0 0 0 0 0]');
+n=size(t,1);
+v=[];d=[];f=[];miss=[];
+for i=1:n
+v(i,:)=norm(x(i,10:12));
+end
+[u,D,fT]=missile_state_u(t,x);
+for i=1:n-1
+d(i,:)=norm(D(i,:)); % drag
+f(i,:)=norm(fT(:,i)); % thrust magnitude
+miss(i,:)=norm(x(i,7:9)-x(i,1:3)); % miss distance
+end

Advanced_control_of_aircraft_spacecraft_Tewari/4.8.m

@@ -0,0 +1,46 @@
+% Calling program for solving the gravity turn 2PBVP
+% by collocation method by MATLAB's intrinsic code 'bvp4c.m'.
+% Requires 'fiveode.m' and 'fivebc.m'.
+% (c) 2009 Ashish Tewari
+% dy/dx=f(y,x); a<=x<=b
+% y(x=a), y(x=b): Boundary conditions
+% y(1,1)=r (km)
+% y(2,1)=v (km/s)
+% y(3,1)= phi (rad)
+% y(4,1)=lambda_2 (Lagrange multiplier of r)
+% y(5,1)=lambda_3 (Lagrange multiplier of v)
+global mu; mu=398600.4; % Grav. const. (km^3/s^2)
+global r0; r0=6378.14; % Planetary radius (km)
+global R;
+R=1; % Control cost coefficient
+tf=600; % Terminal time (s)
+% Collocation points & initial guess follow:
+solinit = bvpinit(linspace(0,tf,5),[6680 8 0.1 0 0]);
+% 2PBVP Solution by collocation method:
+sol = bvp4c(@fiveode,@fivebc,solinit);
+x = linspace(0,tf); % Time vector (s)
+y = deval(sol,x); % Solution state vector
+
+% Program for specifying governing ODEs expressed as
+% state equations for the 2PBVP (to be called by 'bvp4c.m')
+function dydx=fiveode(x,y)
+global mu;
+global r0;
+global R;
+% State equations follow:
+dydx=[y(2)*sin(y(3))
+-0.5*y(5)/R-mu*sin(y(3))/y(1)^2
+(y(2)-mu/(y(1)*y(2)))*cos(y(3))/y(1)
+-2*mu*y(5)*sin(y(3))/y(1)^3
+-y(4)*sin(y(3))];
+% Program for specifying boundary conditions for the 2PBVP.
+% (To be called by 'bvp4c.m')
+function res=fivebc(ya,yb)
+global mu;
+global r0;
+% Boundary conditions follow:
+res=[ya(1)-r0
+ya(2)-0.001
+yb(1)-r0-300
+yb(2)-sqrt(mu/(r0+300))
+yb(3)];

Advanced_control_of_aircraft_spacecraft_Tewari/4.9.m

@@ -0,0 +1,51 @@
+% Calling program for solving the gravity turn,
+% constant fwd acceleration 2PBVP by collocation
+% method using MATLAB's intrinsic code 'bvp4c.m'.
+% Requires 'constuode.m' and 'constubc.m'.
+% (c) 2009 Ashish Tewari
+% dy/dx=f(y,x); a<=x<=b
+% y(x=a), y(x=b): Boundary conditions
+% y(1,1)=r (km)
+% y(2,1)=v (km/s)
+% y(3,1)=phi (rad)
+global mu; mu=398600.4; % Grav. const. (km^3/s^2)
+global r0; r0=6378.14; % Planetary radius (km)
+global R; R=3.6633e6; % Control cost coefficient
+global tf; tf=562; % Terminal time (s)
+% Collocation points & initial guess follow:
+solinit = bvpinit(linspace(0,tf,5),[6680 8 0.1]);
+% 2PBVP Solution by collocation method:
+sol = bvp4c(@constuode,@constubc,solinit);
+x = linspace(0,tf); % Time vector (s)
+y = deval(sol,x); % Solution state vector
+
+% Program for specifying governing state equations
+% for the 2PBVP (to be called by 'bvp4c.m')
+function dydx=constuode(x,y)
+global tf;
+global mu;
+global r0;
+global R;
+if y(1)>=r0
+if x<tf
+u=sqrt(2*tf/R);
+else
+u=0;
+end
+% State equations follow:
+dydx=[y(2)*sin(y(3))
+u-mu*sin(y(3))/y(1)^2
+(y(2)-mu/(y(1)*y(2)))*cos(y(3))/y(1)];
+else
+dydx=[0 0 0]';
+end
+
+% Program for specifying boundary conditions for the 2PBVP.
+% (To be called by 'bvp4c.m')
+function res=constubc(ya,yb)
+global mu;
+global r0;
+% Boundary conditions follow:
+res=[ya(1)-r0
+     ya(2)-0.001
+     yb(3)];

Advanced_control_of_aircraft_spacecraft_Tewari/6.4.m

@@ -0,0 +1,42 @@
+% (c) 2010 Ashish Tewari
+global mu; mu=398600.4e9;
+global R0; R0=6378.14e3;
+global omega; omega = 2*pi/(23*3600+56*60+4.0905);%rad/sec
+global hf; hf=120e3;
+global delf; delf=26*pi/180;
+global lamf; lamf=-82*pi/180;
+global S; S =4; %reference base area, m2
+global c; c=0.5; % nose-radius, m
+global tb1; tb1=100; %first-stage burn-time (s)
+global tb2; tb2=175; %second-stage burn-time (s)
+global fT1;fT1=824532.879800791/2; %first-stage thrust (N)
+global fT2;fT2=161304.620971835*(1.6/1.75); %second-stage thrust (N)
+global m01; m01=30856.5129807023; %first-stage initial mass (kg)
+global m02; m02=8262.36174702614; %second-stage initial mass (kg)
+global mp1; mp1=21012.5606473188; %first-stage propellant mass (kg)
+global mp2; mp2=7516.74365967484; %second-stage propellant mass (kg)
+global mL; mL=350; %payload mass (kg)
+global Gamma; Gamma=1.41; %specific heat ratio, cp/cv
+global machr; machr = [0 0.2 0.3 0.4 0.5 0.6 0.8 0.9 0.95 1.05 1.1 ...
+1.2 1.6 2.0 2.5 3 3.8 5 10 99];
+global Cdr; Cdr =[.475475 .475475 .47576 .48336 .488965 .508345 ...
+.56563 .618165 .668135 1.031795 1.01707 .990565 ...
+.815955 .69236 .60971 .54606 .513 .494 .48317 .48317];
+dtr=pi/180;
+%Initial conditions:
+long = -80.55*dtr;
+%initial longitude.
+lat = 28.5*dtr;
+%initial latitude.
+rad=R0;
+%initial radius (m)
+vel=0.0001;
+%initial velocity (m/s)
+fpa=90*dtr;
+%initial f.p.a
+chi=170*dtr; %initial heading angle (measured from north)
+orbinit = [long; lat; rad; vel; fpa; chi];
+[t,x]=ode45('lambert_rocket',[0 275],orbinit);
+er=(x(size(x,1),3)-6378.14e3-hf)/1000 %final radial error (km)
+edel=(x(size(x,1),2)-delf)/dtr %final latitude error (deg.)
+elam=(x(size(x,1),1)-lamf)/dtr %final longitude error (deg.)

Advanced_control_of_aircraft_spacecraft_Tewari/airnav.m

@@ -0,0 +1,13 @@
+% Program for specifying governing ODEs of great circle navigation
+% expressed as state equations
+% (c) 2010 Ashish Tewari
+function dydx=airnav(x,y)
+global vprime;
+global vw;
+global Aw;
+global r;
+global Omega;
+v=sqrt(vprime^2-vw^2*sin(Aw-y(3))^2)+vw*cos(Aw-y(3));
+dydx=[v*cos(y(3))/r
+v*sin(y(3))/(r*cos(y(1)))
+v*sin(y(3))*tan(y(1))/r];

Advanced_control_of_aircraft_spacecraft_Tewari/airnavode.m

@@ -0,0 +1,41 @@
+% Program for specifying governing ODEs for optimal navigation
+% expressed as state equations for the 2PBVP
+%(to be called by 'bvp4c.m')
+% (c) 2010 Ashish Tewari
+function dydx=airnavode(x,y)
+global vprime;
+global vw;
+global Aw;
+global r;
+global R;
+global Omega;
+v=sqrt(vprime^2-vw^2*sin(Aw-y(3))^2)+vw*cos(Aw-y(3));
+u=-0.5*y(6)/(R*v);
+dvdA=sin(Aw-y(3))*(1-cos(Aw-y(3))/sqrt(vprime^2-vw^2*sin(Aw-y(3))^2));
+dydx=[v*cos(y(3))/r
+v*sin(y(3))/(r*cos(y(1)))
+v*sin(y(3))*tan(y(1))/r+u/v+Omega^2*r*sin(y(3))*sin(2*y(1))/(2*v)+ ...
+2*Omega*sin(y(1))
+v*sin(y(3))*tan(y(1))*sec(y(1))*y(5)/r+(v*sin(y(3))*sec(y(1))^2/r+ ...
+Omega^2*r*sin(y(3))*cos(2*y(1))/v+2*Omega*cos(y(1)))*y(6)
+0
+(dvdA*cos(y(3))-v*sin(y(3)))*y(4)/r + ...
+(v*cos(y(3))+dvdA*sin(y(3)))*y(5)/(r*cos(y(1))) + ...
+(tan(y(1))*(v*cos(y(3))+dvdA*sin(y(3)))/r-u*dvdA/v^2+ ...
+ 0.5*Omega^2*r*sin(2*y(1))*(cos(y(3))/v-dvdA*sin(y(3))/v^2))*y(6)];
+
+% Program for specifying boundary conditions for 2PBVP of
+% optimal aircraft navigation.
+% (To be called by bvp4c.m)
+function res=airnavbc(ya,yb)
+global lam0;
+global lamf;
+global del0;
+global delf;
+global A0;
+res=[ya(1)-del0
+ya(2)-lam0
+ya(3)-A0
+yb(1)-delf
+yb(2)-lamf
+yb(6)];