攻击时间控制的动态逆三维制导律源代码
2019-04-13 21:47发布
生成海报
%参考文献:《攻击时间控制的动态逆三维制导律》,哈尔滨工程大学学报,2010
%编程难点:视线坐标系与参考坐标系间的转换关系
clear
clc
V=300;% 导弹速度,m/s
x=-5000;%进攻弹坐标,m
y=5000;
z=6000;
xt=0;%目标坐标,m
yt=0;
zt=0;
Rm=sqrt((x-xt)^2+(y-yt)^2+(z-zt)^2);
DRm=0;
theta_m=20*pi/180;
Dtheta_m=0;
fea_m=10*pi/180;
theta_s=atan((zt-z)/sqrt((xt-x)^2+(yt-y)^2));
fea_s=atan((yt-y)/(xt-x));
Dtheta_s=0;
Dfea_s=0;
g=9.8;%重力加速度
k1=10;
c1=0.7;
c2=0.9;
k3=5;
Td=45;
n=1;
t=0;
dt=0.001;
while DRm<=0 &&Rm>=0
et=Td-t-Rm/V;
k2=(1-cos(theta_m)+c1*et*cos(theta_m)/(et+c2))/et;
Aym=k1*abs(DRm)*Dtheta_s/cos(theta_m);%增广比例导引律
if Aym>8*g
Aym=8*g;
end
if Aym<-8*g
Aym=-8*g;
end
if fea_m>=0 && fea_m=-pi/2 && fea_m<0
Azm=-(V*(1-k2*et)*sin(theta_m)*Dtheta_m/cos(theta_m)+k2*V-k2*V*cos(theta_m)*cos(fea_m))/sqrt(1-((1-k2*et)/cos(theta_m))^2)+k3*V*cos(theta_m)*(fea_m+acos((1-k2*et)/cos(theta_m)))+V^2/Rm*sin(fea_m)*(1-sin(theta_m)*cos(theta_m)*cos(fea_m)*tan(theta_s));
end
% Azm=-k1*abs(DRm)*Dfea_s/cos(fea_m);%方法2:偏航通道也采用增广比例导引律
if Azm>8*g
Azm=8*g;
end
if Azm<-8*g
Azm=-8*g;
end
DRm=-V*cos(theta_m)*cos(fea_m);
Dtheta_m=Aym/V+V*cos(theta_m)*sin(fea_m)*sin(fea_m)*tan(theta_s)/Rm+V*sin(theta_m)*cos(fea_m)/Rm;
Dfea_m=-Azm/(V*cos(theta_m))+V*sin(fea_m)/(Rm*cos(theta_m))-V*sin(theta_m)*cos(fea_m)*tan(theta_s)*sin(fea_m)/Rm;
Dtheta_s=-V*sin(theta_m)/Rm;
Dfea_s=-V*cos(theta_m)*sin(fea_m)/(Rm*cos(theta_s));
Rm=Rm+DRm*dt;
theta_m=theta_m+Dtheta_m*dt;
fea_m=fea_m+Dfea_m*dt;
theta_s=theta_s+Dtheta_s*dt;
fea_s=fea_s+Dfea_s*dt;
%视线坐标系转换到参考坐标系-
Tran=[cos(theta_s)*cos(fea_s) sin(theta_s) -cos(theta_s)*sin(fea_s);
-sin(theta_s)*cos(fea_s) cos(theta_s) sin(theta_s)*sin(fea_s);
sin(fea_s) 0 cos(fea_s)];
DP=Tran'*[V*cos(theta_m)*cos(fea_m);V*cos(theta_m)*sin(fea_m);V*sin(theta_m)];
Dx=DP(1);
Dy=DP(2);
Dz=DP(3);
x=x+Dx*dt;
y=y+Dy*dt;
z=z+Dz*dt;
P_s(:,n)=[x;y;z];
S_s(:,n)=[Rm;theta_m;fea_m;et];
A_s(:,n)=[Aym/g;Azm/g];
n=n+1;
t=t+dt;
end
figure(1)%由于plot3函数默认采用美国坐标系,需要转化为我们习惯的俄式坐标系
plot3(P_s(2,:)/1000,P_s(1,:)/1000,-P_s(3,:)/1000);
xlabel('Y/km')
ylabel('X/km')
zlabel('Z/km')
figure(2)
plot((1:n-1)*dt,S_s(1,:))
xlabel('time/s')
ylabel('相对距离')
figure(3)
plot((1:n-1)*dt,S_s(2,:)*180/pi)
xlabel('time/s')
ylabel('俯仰通道前置角/deg')
figure(4)
plot((1:n-1)*dt,S_s(3,:)*180/pi)
xlabel('time/s')
ylabel('偏航通道前置角/deg')
figure(5)
plot((1:n-1)*dt,S_s(4,:))
xlabel('time/s')
ylabel('et/s')
figure(6)
plot((1:n-1)*dt,A_s(1,:))
xlabel('time/s')
ylabel('Amy/g')
title('俯仰通道控制量')
figure(7)
plot((1:n-1)*dt,A_s(2,:))
xlabel('time/s')
ylabel('Amz/g')
title('偏航通道控制量')
打开微信“扫一扫”,打开网页后点击屏幕右上角分享按钮