全国大学生数学建模竞赛2014A题嫦娥三号软着陆轨道设计与控制策略MATLAB程序
程序员文章站
2022-03-31 23:26:03
...
目录
1.3 不考虑飞船飞行角度调整时计算轨迹并绘图MATLAB程序
2.1 考虑飞船飞行角度调整时计算轨迹并绘图MATLAB程序
一、第1问
1.1 图像预处理MATLAB程序
%% 图像处理 A=imread('附件3 距2400m处的数字高程图.tif'); surf(A); set(gca,'FontSize',22); shading interp; colorbar; colormap gray; figure B=imread('附件4 距月面100m处的数字高程图.tif'); mesh(B); set(gca,'FontSize',22); shading interp; colorbar; colormap gray;
1.2 龙格库塔法求解微分方程初始化函数MATLAB程序
%微分方程组 function df=odefun(t,y,ai) G=6.67259e-11;%引力常量 M=7.3477e22;%月球质量 m0=2400;%飞行器质量 r=1737.013e3;%月球半径 a=15e3;%近月点距离 b=100e3;%远月点距离 Fmax=7500;%最大推力 rA=r+a;%椭圆短半轴 rB=r+b;%椭圆长半轴 VB=sqrt(2*G*M*(rA/(rB*(rA+rB))));%远月点速度 VA=(rB/rA)*VB;%近月点速度 df=zeros(4,1); yy=y(1);dy=y(2);xx=y(3);dx=y(4); m=7500/3087*t; df(1)=dy; df(2)=G*M/(yy^2)-Fmax*sind(ai)/(m0-m); df(3)=dx; df(4)=-Fmax*cosd(ai)/(m0-m); % df(1)=dy; % df(2)=-G*M/(xx^2+yy^2)*cosd(atand(xx/yy))+Fmax*sind(ai)/(m0-m); % df(3)=dx; % df(4)=-G*M/(xx^2+yy^2)*sind(atand(xx/yy))-Fmax*cosd(ai)/(m0-m); end
1.3 不考虑飞船飞行角度调整时计算轨迹并绘图MATLAB程序
%% 求速度 G=6.67259e-11;%引力常量 M=7.3477e22;%月球质量 r=1737.013e3;%月球半径 a=15e3;%近月点距离 b=100e3;%远月点距离 R=r+a-2640-12000;%最终海拔 V=57;%最终垂直速度 rA=r+a;%椭圆短半轴 rB=r+b;%椭圆长半轴 VB=sqrt(2*G*M*(rA/(rB*(rA+rB))));%远月点速度 VA=(rB/rA)*VB;%近月点速度 %% 解微分方程（遍历ai） tspan=[0 500]; y0=[rA 0 0 VA]; ri=zeros(93,1);%寻找最接近海拔R的高度 vi=zeros(93,1);%寻找最接近速度V的值 m1t=zeros(11,1);%寻找满足R的时间的值 m2t=zeros(11,1);%寻找满足V的时间的值 count=0;%循环计数 for ai=25:35 count=count+1; [t,y]=ode45(@(t,y)odefun(t,y,ai),tspan,y0); for i=1:93 ri(i)=abs(y(i,1)-R); vi(i)=abs(sqrt((y(i,2))^2+(y(i,4))^2)-V); end ii1=find(ri==min(ri)); ii2=find(vi==min(vi)); m1t(count)=t(ii1(1));%20个ai角度对应满足R的时间的值 m2t(count)=t(ii2(1));%20个ai角度对应满足V的时间的值 end %% 解微分方程（代入ai求解画图） tspan=[0 500]; y0=[rA 0 0 VA]; ri=zeros(93,1);%寻找最接近海拔R的高度 vi=zeros(93,1);%寻找最接近速度V的值 m1t=zeros(11,1);%寻找满足R的时间的值 m2t=zeros(11,1);%寻找满足V的时间的值 count=0;%循环计数 for ai=29 %遍历得到最优解29° count=count+1; [t,y]=ode45(@(t,y)odefun(t,y,ai),tspan,y0); for i=1:93 ri(i)=abs(y(i,1)-R); vi(i)=abs(sqrt((y(i,2))^2+(y(i,4))^2)-V); end ii1=find(ri==min(ri)); ii2=find(vi==min(vi)); m1t(count)=t(ii1(1));%20个ai角度对应满足R的时间的值 m2t(count)=t(ii2(1));%20个ai角度对应满足V的时间的值 end L=0;%总长度 for ii=1:ii1 L=L+sqrt((y(ii,2))^2+(y(ii,4))^2)*(t(ii+1)-t(ii)); end si=44.12-L/r*180/pi;%纬度（经度不变） plot(t,y(:,1),'linewidth',2) axis([0 450 1.725e6 1.755e6]); xlabel('t(s)','FontSize',28); ylabel('y(m)','FontSize',28); set(gca,'FontSize',28,'linewidth',1);
二、第2问
2.1 考虑飞船飞行角度调整时计算轨迹并绘图MATLAB程序
%% 求速度 G=6.67259e-11;%引力常量 M=7.3477e22;%月球质量 r=1737.013e3;%月球半径 a=15e3;%近月点距离 b=100e3;%远月点距离 R=r+a-2640-12000;%最终海拔 V=57;%最终垂直速度 rA=r+a;%椭圆短半轴 rB=r+b;%椭圆长半轴 VB=sqrt(2*G*M*(rA/(rB*(rA+rB))));%远月点速度 VA=(rB/rA)*VB;%近月点速度 %% 解微分方程（第一次遍历ai） tt=zeros(1,1);%记录满足条件的t5值 flag=0;%满足条件计数 tspan1=[0 100]; y01=[rA 0 0 VA]; count1=0;%循环计数 for ai1=20:40 count1=count1+1; [t1,y1]=ode45(@(t,y)odefun(t,y,ai1),tspan1,y01); n1=size(t1,1); %% 解微分方程（第二次遍历ai） tspan2=[100 200]; y02=[y1(n1,1) y1(n1,2) y1(n1,3) y1(n1,4)]; count2=0;%循环计数 for ai2=20:40 count2=count2+1; [t2,y2]=ode45(@(t,y)odefun(t,y,ai2),tspan2,y02); n2=size(t2,1); %% 解微分方程（第三次遍历ai） tspan3=[200 300]; y03=[y2(n2,1) y2(n2,2) y2(n2,3) y2(n2,4)]; count3=0;%循环计数 for ai3=10:30 count3=count3+1; [t3,y3]=ode45(@(t,y)odefun(t,y,ai3),tspan3,y03); n3=size(t3,1); %% 解微分方程（第四次遍历ai） tspan4=[300 400]; y04=[y3(n3,1) y3(n3,2) y3(n3,3) y3(n3,4)]; count4=0;%循环计数 for ai4=10:30 count4=count4+1; [t4,y4]=ode45(@(t,y)odefun(t,y,ai4),tspan4,y04); n4=size(t4,1); %% 解微分方程（第五次遍历ai） tspan5=[400 500]; y05=[y4(n4,1) y4(n4,2) y4(n4,3) y4(n4,4)]; count5=0;%循环计数 for ai5=20:40 count5=count5+1; [t5,y5]=ode45(@(t,y)odefun(t,y,ai5),tspan5,y05); n5=size(t5,1); flag=flag+1; for i=1:n5 if (y5(i,1)-R<0)&&(sqrt((y5(i,2))^2+(y5(i,4))^2)-V<0) tt(flag,1)=t5(i); break end end end end end end end %% 数据处理 ind1=find(tt==0);%找出tt中所有为0的标号 tt(ind1)=inf; ind2=find(tt==400);%找出tt中所有为400的标号 tt(ind2)=inf; ind0=find(tt==min(tt));%找出tt中最小值 mi=tt(ind0); ind0=ind0(1); mi=mi(1); %% 遍历求得最优最短时间后代入 %% 求速度 G=6.67259e-11;%引力常量 M=7.3477e22;%月球质量 r=1737.013e3;%月球半径 a=15e3;%近月点距离 b=100e3;%远月点距离 R=r+a-2640-12000;%最终海拔 V=57;%最终垂直速度 rA=r+a;%椭圆短半轴 rB=r+b;%椭圆长半轴 VB=sqrt(2*G*M*(rA/(rB*(rA+rB))));%远月点速度 VA=(rB/rA)*VB;%近月点速度 %% 解微分方程（第一次遍历ai） tt=zeros(1,1);%记录满足条件的t5值 flag=0;%满足条件计数 tspan1=[0 100]; y01=[rA 0 0 VA]; count1=0;%循环计数 for ai1=20:40 count1=count1+1; [t1,y1]=ode45(@(t,y)odefun(t,y,ai1),tspan1,y01); n1=size(t1,1); %% 解微分方程（第二次遍历ai） tspan2=[100 200]; y02=[y1(n1,1) y1(n1,2) y1(n1,3) y1(n1,4)]; count2=0;%循环计数 for ai2=20:40 count2=count2+1; [t2,y2]=ode45(@(t,y)odefun(t,y,ai2),tspan2,y02); n2=size(t2,1); %% 解微分方程（第三次遍历ai） tspan3=[200 300]; y03=[y2(n2,1) y2(n2,2) y2(n2,3) y2(n2,4)]; count3=0;%循环计数 for ai3=10:30 count3=count3+1; [t3,y3]=ode45(@(t,y)odefun(t,y,ai3),tspan3,y03); n3=size(t3,1); %% 解微分方程（第四次遍历ai） tspan4=[300 400]; y04=[y3(n3,1) y3(n3,2) y3(n3,3) y3(n3,4)]; count4=0;%循环计数 for ai4=10:30 count4=count4+1; [t4,y4]=ode45(@(t,y)odefun(t,y,ai4),tspan4,y04); n4=size(t4,1); %% 解微分方程（第五次遍历ai） tspan5=[400 500]; y05=[y4(n4,1) y4(n4,2) y4(n4,3) y4(n4,4)]; count5=0;%循环计数 for ai5=20:40 count5=count5+1; [t5,y5]=ode45(@(t,y)odefun(t,y,ai5),tspan5,y05); n5=size(t5,1); flag=flag+1; if flag==ind0 break end end if flag==ind0 break end end if flag==ind0 break end end if flag==ind0 break end end if flag==ind0 break end end %% plot(t1,y1(:,1),'b','linewidth',2) axis([0 350 1.725e6 1.755e6]) xlabel('t(s)','FontSize',28); ylabel('y(m)','FontSize',28); set(gca,'FontSize',28,'linewidth',1); hold on plot(t2,y2(:,1),'b','linewidth',2) axis([0 350 1.725e6 1.755e6]) xlabel('t(s)','FontSize',28); ylabel('y(m)','FontSize',28); set(gca,'FontSize',28,'linewidth',1); hold on plot(t3,y3(:,1),'b','linewidth',2) axis([0 350 1.725e6 1.755e6]) xlabel('t(s)','FontSize',28); ylabel('y(m)','FontSize',28); set(gca,'FontSize',28,'linewidth',1); hold on plot(t4,y4(:,1),'b','linewidth',2) axis([0 350 1.725e6 1.755e6]) xlabel('t(s)','FontSize',28); ylabel('y(m)','FontSize',28); set(gca,'FontSize',28,'linewidth',1); hold on plot(t5,y5(:,1),'b','linewidth',2) axis([0 350 1.725e6 1.755e6]) xlabel('t(s)','FontSize',28); ylabel('y(m)','FontSize',28); set(gca,'FontSize',28,'linewidth',1);
2.2 实现粗避障MATLAB程序
%% 绘图（粗避障） A=imread('附件3 距2400m处的数字高程图.tif'); A=double(A); n=size(A,1); for i=1:n for j=1:n if A(i,j)<80||A(i,j)>100 A(i,j)=0; end end end surf(A) set(gca,'FontSize',22); shading interp; colorbar; colormap gray; %% 计算方差、距离 V=zeros(n-100,n-100);%方差 D=zeros(n-100,n-100);%距离 for i=1:n-100 for j=1:n-100 AA=A(i:99+i,j:99+j); if find(AA==0) V(i,j)=inf; D(i,j)=inf; else V(i,j)=var(AA(:)); D(i,j)=sqrt( (((99+i+i)/2)-n/2)^2+(((99+j+j)/2)-n/2)^2 ); end end end count=size(find(V==inf),1);%找出不符合着陆要求的点的个数 count=(n-1)^2-count;%计算符合要求的点的个数 sum1=zeros(count,1);%定义用于求方差和 sum2=zeros(count,1);%定义用于求距离和 count0=0;%用于循环计数 for i=1:n-100 for j=1:n-100 if V(i,j)==inf||D(i,j)==inf continue else count0=count0+1; sum1(count0)=V(i,j); sum2(count0)=D(i,j); end end end V0=mean(sum1);%平均方差 D0=mean(sum2);%平均距离 mi=zeros(n-100,n-100);%加权方差与距离 for i=1:n-100 for j=1:n-100 if V(i,j)==inf||D(i,j)==inf continue else mi(i,j)=V(i,j)/V0+D(i,j)/D0; end end end %% 找最小值 ind=find(mi==0);%找出mi中所有为0的位置标号 mi(ind)=inf;%用无穷代替 [x,y]=find(mi==min(min(mi))); m=mi(x,y); %% 画方框图 hold on plot3([y y+99],[x x],[100 100],'r','linewidth',2) hold on plot3([y y],[x x+99],[100 100],'r','linewidth',2) hold on plot3([y y+99],[x+99 x+99],[100 100],'r','linewidth',2) hold on plot3([y+99 y+99],[x x+99],[100 100],'r','linewidth',2)
2.3 实现细避障MATLAB程序
%% 绘图（细避障） B=imread('附件4 距月面100m处的数字高程图.tif'); B=double(B); n=size(B,1); for i=1:n for j=1:n if B(i,j)<80||B(i,j)>130 B(i,j)=0; end end end surf(B) set(gca,'FontSize',22); shading interp; colorbar; colormap gray; %% 计算方差、距离 V=zeros(n-30,n-30);%方差 D=zeros(n-30,n-30);%距离 for i=1:n-30 for j=1:n-30 BB=B(i:29+i,j:29+j); if find(BB==0) V(i,j)=inf; D(i,j)=inf; else V(i,j)=var(BB(:)); D(i,j)=sqrt( (((29+i+i)/2)-n/2)^2+(((29+j+j)/2)-n/2)^2 ); end end end count=size(find(V==inf),1);%找出不符合着陆要求的点的个数 count=(n-1)^2-count;%计算符合要求的点的个数 sum1=zeros(count,1);%定义用于求方差和 sum2=zeros(count,1);%定义用于求距离和 count0=0;%用于循环计数 for i=1:n-30 for j=1:n-30 if V(i,j)==inf||D(i,j)==inf continue else count0=count0+1; sum1(count0)=V(i,j); sum2(count0)=D(i,j); end end end V0=mean(sum1);%平均方差 D0=mean(sum2);%平均距离 mi=zeros(n-30,n-30);%加权方差与距离 for i=1:n-30 for j=1:n-30 if V(i,j)==inf||D(i,j)==inf continue else mi(i,j)=V(i,j)/V0+D(i,j)/D0; end end end %% 找最小值 ind=find(mi==0);%找出mi中所有为0的位置标号 mi(ind)=inf;%用无穷代替 [x,y]=find(mi==min(min(mi))); m=mi(x,y); %% 画方框图 hold on plot3([y y+29],[x x],[130 130],'r','linewidth',2) hold on plot3([y y],[x x+29],[130 130],'r','linewidth',2) hold on plot3([y y+29],[x+29 x+29],[130 130],'r','linewidth',2) hold on plot3([y+29 y+29],[x x+29],[130 130],'r','linewidth',2)