高超声速飞行器俯仰通道滑模控制律设计。
代码如下：
clc;clear all;%********模型参数******&*%gamma=0;m=5640*14.593903;%质量Iz=6.07*14.593903/(0.3048)^2*10^6;%转动惯量s=2690*(0.3048)^2;%参考面积L=39.6*0.3048;%特征长度%模型初始条件V=3400;%速度Ma=5;h=36000;%高度Re=6371004;%地球半径mu=3.986005*10^14;%万有引力常数r=h+Re;%地心距g=mu/r/r;rho=0.36392/exp((h-11000)/6341.62);%大气密度q=0.5*rho*V*V;%动压alpha=[];%攻角delta=[];%舵偏角wz=[];%角速度t=[];%时间dt=0.1;%步长%气动系数Cy=[];Calpha=[];Cdelta=[];CY=[];Cmalpha=[];Cmdelta=[];Cmwz=[];CM=[];Y=[];%升力Mz=[];%力矩%*********控制器参数***********%c=20;k=1;epsilon=1;nyc=[];%过载指令dnyc=[];%过载指令变化率ddnyc=[];%过载指令二阶ny=[];%实际过载dny=[];%实际过载变化率e=[];de=[];S=[];a1=[];a2=[];deltac=[];%*******初始值**********%%*********模型部分************%alpha(1)=0;delta(1)=0;wz(1)=0;%*******气动参数*************%alpha0=[12.5*pi/180,16*pi/180,20.5*pi/180,25*pi/180,29*pi/180,35*pi/180,38*pi/180,40*pi/180];Y0=[0.4,0.38,0.5,0.57,0.7,0.8,0.87,0.86];Yalpha0=[1.8,1.6,1.63,1.6,1.6,1.4,.0,1.19];PY= polyfit(alpha0,Y0,3);Cy(1)=PY(1)*(alpha(1))^3+PY(2)*(alpha(1))^2+PY(3)*(alpha(1))+PY(4);Palpha= polyfit(alpha0,Yalpha0,3);Calpha(1)=Palpha(1)*(alpha(1))^3+Palpha(2)*(alpha(1))^2+Palpha(3)*(alpha(1))+Palpha(4);CY(1)=Calpha(1)*alpha(1);Y(1)=CY(1)*q*s;Malpha=[-0.16,-0.05,-0.05,-0.075,-0.1,-0.14,-0.15,-0.175];Mdelta=[-0.17,-0.1,-0.11,-0.12,-0.135,-0.16,-0.22,-0.27];Mwz=[-1.1,-0.9,-1.1,-1.1,-1.2,-1.3,-1.5,-1.6];Pmalpha= polyfit(alpha0,Malpha,3);Cmalpha(1)=Pmalpha(1)*(alpha(1))^3+Pmalpha(2)*(alpha(1))^2+Pmalpha(3)*(alpha(1))+Pmalpha(4);Pmdelta= polyfit(alpha0,Mdelta,3);Cmdelta(1)=Pmdelta(1)*(alpha(1))^3+Pmdelta(2)*(alpha(1))^2+Pmdelta(3)*(alpha(1))+Pmdelta(4);Pwz= polyfit(alpha0,Mwz,3);Cmwz(1)=Pwz(1)*(alpha(1))^3+Pwz(2)*(alpha(1))^2+Pwz(3)*(alpha(1))+Pwz(4);CM(1)=Cmdelta(1)*delta(1);Mz(1)=CM(1)*q*s*L;%********控制器初始值*******%a1(1)=Calpha(1)*q*s/(m*V);a2(1)=Cmdelta(1)*q*s*L/Iz;ny(1)=Y(1)/(m*g);dny(1)=a1(1)*wz(1)*V/g-a1(1)*ny(1);ny(1)=Y(1)/(m*g);nyc(1)=ny(1);dnyc(1)=0;ddnyc(1)=0;e(1)=nyc(1)-ny(1);de(1)=dnyc(1)-dny(1);S(1)=c*e(1)+de(1);delta(1)=0;%************仿真**************%t(1)=0;i=2;while (i<=1000)%*************控制器模型************%a1(i)=Calpha(i-1)*q*s/(m*V);a2(i)=Cmdelta(i-1)*q*s*L/Iz;nyc(i)=sin(i*dt*0.1);dnyc(i)=(nyc(i)-nyc(i-1))/dt;ddnyc(i)=(dnyc(i)-dnyc(i-1))/dt;ny(i)=Y(i-1)/(m*g);dny(i)=a1(i-1)*wz(i-1)*V/g-a1(i-1)*ny(i-1);e(i)=nyc(i)-ny(i-1);de(i)=dnyc(i)-dny(i-1);S(i)=c*e(i)+de(i);delta(i)=(k*S(i)+epsilon*sign(S(i))+c*dnyc(i)+ddnyc(i)+V/g*a1(i)*...(a1(i)-c)*wz(i-1)+a1(i)*(c-a1(i))*ny(i))/(V*a1(i)*a2(i)/g);if (delta(i)>25/180*pi)delta(i)=25/180*pi;else if(delta(i)<-25/180*pi)delta(i)=-25/180*pi;endend%*********俯仰通道模型*********%alpha(i)=(wz(i-1)-Y(i-1)/(m*V)-(mu-V*V*r)*cos(gamma)/(V*r*r))*dt+alpha(i-1);wz(i)=Mz(i-1)/Iz*dt+wz(i-1);Calpha(i)=Palpha(1)*(alpha(i))^3+Palpha(2)*(alpha(i))^2+Palpha(3)*(alpha(i))+Palpha(4);CY(i)=Calpha(i)*alpha(i);Y(i)=CY(i)*q*s;Cmalpha(i)=Pmalpha(1)*(alpha(i))^3+Pmalpha(2)*(alpha(i))^2+Pmalpha(3)*(alpha(i))+Pmalpha(4);Cmdelta(i)=Pmdelta(1)*(alpha(i))^3+Pmdelta(2)*(alpha(i))^2+Pmdelta(3)*(alpha(i))+Pmdelta(4);Cmwz(i)=Pwz(1)*(alpha(i))^3+Pwz(2)*(alpha(i))^2+Pwz(3)*(alpha(i))+Pwz(4);CM(i)=Cmalpha(i)*alpha(i)+Cmdelta(i)*delta(i)+Cmwz(i)*wz(i);Mz(i)=CM(i)*q*s*L;t(i)=t(i-1)+dt;i=i+1;endfigure(1);plot(t,ny,'*');hold on;xlabel("时间/s");ylabel("过载");plot(t,nyc);legend("飞行器实际过载","过载指令");
其中气动参数部分有误。
