2014年国赛高教杯数学建模
A题 嫦娥三号软着陆轨道设计与控制策略
嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m（见附件1）。
嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段（见附件2），要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。
根据上述的基本要求，请你们建立数学模型解决下面的问题：
（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。
（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
（3）对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
附件1： 问题的背景与参考资料；
附件2： 嫦娥三号着陆过程的六个阶段及其状态要求；
附件3：距月面2400m处的数字高程图；
附件4：距月面100m处的数字高程图。
整体求解过程概述(摘要)
本文研究的是嫦娥三号软着陆轨道设计与最优控制策略问题，根据动力学相关原理，建立了嫦娥三号软着陆轨迹模型，以燃料消耗量最小为目标，得到软着陆过程中各阶段的最优控制策略。
针对问题1，首先确定通过嫦娥三号着陆准备轨道位于其着陆点所在经线与月心构成的平面内，然后对嫦娥三号在主减速段进行受力分析，以近月点在月球表面的投影点为原点，以着陆轨道参考系纵向平面为xoy面，以初始飞行方向为x轴，以原点与近月点的连线方向为y轴，建立平面直角坐标系，根据牛顿第二定律，建立嫦娥三号主减速段轨迹模型，以燃料消耗量最小为目标，通过轨迹离散化，逐步迭代求得该阶段的水平位移，再依据地理学经纬度计算规则，建立地表距离与经纬度转化模型，最终得到近月点在月球表面的投影位置为19.51’W,31.68’N，距离月球表面高度为15km，远月点在月球表面的投影位置为160.49’E,31.68’S，距离月球表面高度为100km。利用机械能守恒定理和开普勒第二定律，最终得到近月点与远月点速度的大小分别为为1.692 km/s ,1.614km/s。
针对问题2，根据牛顿第二定律，以每个阶段初始点和终值点的状态为约束，以燃料消耗最小为目标，建立全局最优模型，通过轨迹离散化，逐步迭代求得每个阶段的水平位移，分别得到软着陆过程6个阶段着陆轨迹方程及其对应的最优控制策略；在粗避障段和精避障段，我们将所给数字高程图均分成9块，综合相对高程差与标准差定义平坦度指标来衡量每一块区域，从而选取最佳着陆点；在粗避障段，分别从两种运动状态：其一是先把主减速发动机关闭，在进行一段时间匀加速直线运动后，再次打开发动机，进行减速直线运动，其二是整个阶段都进行匀速直线运动，以燃料消耗最小为目标，确定前者更优从而确定该阶段最优控制策略。
针对问题3，从改变近月点离月球表面的距离和主减速发动机提供的推力两个方面，对嫦娥三号在该段的水平位移、燃料消耗等参数进行灵敏度分析，发现近月点离月球表面的距离与该段的水平位移和燃料消耗呈线性正相关，发现主减速发动机提供的推力与该段的水平位移呈线性负相关，与该段的燃料消耗呈线性正相关。由于嫦娥三号在主减速段水平位移最大，选取该段从对近月点离月球表面的距离变化和主减速发动机提供的推力变化两个角度对模型进行局部阶段误差分析，通过计算每个阶段时间的相对误差对模型进行整体误差分析。
最后，对模型的优缺点进行评价，并提出改进的方向。
模型假设：
（1）假设月球引力场为垂直于月面的均匀引力场；
（2）假设在几百秒范围内的下降时间里，月球引力非球项、日月引力摄动和月球自转均可以忽略；
（3）假设只考虑惯性和引力作用下的稳定的椭圆运动状态，而对于“嫦娥三号”的变轨状态，情况过于复杂，不予以考虑；
（4）制动发动机的最大推力与初始质量之比大于月面引力加速度,并且制动推进系统能够在一定的初始条件下将探测器停于月面上。
问题分析：
嫦娥三号从实施近月制动到成功着陆主要经历了环月轨道→椭圆轨道→着陆轨道三个变轨过程，从环月轨道下降到着陆点的过程，称为软着陆过程，其包括着陆准备轨道、主减速段、快速调整段、粗避障段、精避障段、缓速下降阶段六个阶段。软着陆：指嫦娥三号以相对月球较小的速度着陆，速度一般为几米/每秒； 比冲：火箭发动机单位质量推进剂产生的冲量，或单位流量的推进剂产生的推力。
针对问题1，首先对软着陆过程进行简化分析，可知主减速段终点已基本位于着陆点上方，其空间坐标在月面上投影即为着陆点坐标，仅考虑主减速段着陆过程，可利用主减速段终点逆推出近月点位置。因此，可取嫦娥三号处于主减速段终点时在月球表面的投影点作为原点，软着陆运动轨迹所在平面建立二维坐标系。结合动力学知识，建立最优主减速轨迹模型，并需满足主减速段燃料消耗最小，得到其着陆轨迹微分方程，利用仿真求解得到近月点坐标，运用地表距离与经纬度转化关系，最终得到近月点与远月点地理位置。结合开普勒第二定律和机械能守恒定律，建立两个相关等式，求得近月点、远月点处速度大小，基于最优主减速段轨迹模型求解结果，便可根据近月点运动方向与坐标轴夹角得到速度方向。
针对问题2，确定嫦娥三号着陆轨道及6个阶段的最优控制策略时，始终要满足燃料消耗最小原则。在问题1中已经对近月点及对主减速段的运动情况进行求解，近月点和主减速段终值点的位置、速度及发动机推力大小均已知，在此基础上，给定准备轨道、主减速段最优控制策略。快速调整段主要是对探测器姿态进行调整，采取与主减速同样的建模方法，得到该段质心动力学方程，在满足约束条件及阶段要求下给出具体最优控制策略。对于粗避障段，首先对其数字高程图进行划分，对每个区域的平坦程度进行分析，取最平坦区域作为着陆大致范围。同样需建立动力学模型对运动轨迹进行描述，考虑到可能存在匀速直线运动和先匀加速后又在恒定推力下减速至0两种情况，在可行情况下，选择燃料消耗最少的方案作为最优策略。
针对问题3， 从改变近月点离月球表面的距离和主减速发动机提供的推力两个方面，对嫦娥三号在该段的水平位移、燃料消耗等参数进行灵敏度分析，进而求得近月点离月球表面的距离与该段的水平位移和燃料消耗之间的关系、主减速发动机提供的推力与该段的水平位移及该段的燃料消耗之间的关系；由于嫦娥三号在主减速段水平位移最大，选取该段从对近月点离月球表面的距离变化和主减速发动机提供的推力变化两个角度对模型进行局部阶段误差分析，通过计算每个阶段时间的相对误差对模型进行整体误差分析。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：
%% 整个软着陆过程的运动轨迹
clc;clear;close all;
g=1.633; %月球重力加速度
m0=2.4*10^3;%卫星初始质量
Ve=2940; %比冲
theta=9.654*pi/180;%初速度与水平方向的夹角
F=7500; %推力
V0=1692.464; %近日点初速度
t=0; %初始时间
T=0.1; %时间步长
Vx0=V0*cos(-theta); %水平初速度
Vy0=V0*sin(-theta); %竖直初速度
Ay0=g-F*sin(-theta)/(m0-F/Ve*t);%竖直初加速度
Ax0=-F*cos(-theta)/(m0-F/Ve*t);%水平初加速度
count=0;
X_res=Vx0*t+0.5*Ax0*t^2;
Y_res=Vy0*t+0.5*Ax0*t^2;
Result=[];
h=12000;%主减速阶段的下落距离
%% 迭代求 分解速度和分解位移
while (Y_res<h )
count=count+1;
Vx=Vx0+Ax0*T;
Vy=Vy0+Ay0*T;
V_res=sqrt(Vx^2+Vy^2);
Vx0=Vx;
Vy0=Vy;
X=Vx0*T+0.5*Ax0*T^2;
Y=Vy0*T+0.5*Ay0*T^2;
X_res=X_res+X;
Y_res=Y_res+Y;Time=count*T;Result=[Result;X_res,Y_res,V_res,Time];
SIN=Vy/sqrt(Vy^2+Vx^2);
COS=Vx/sqrt(Vy^2+Vx^2);
Ay=g-F*SIN/(m0-F/Ve*Time);
Ax=-F*COS/(m0-F/Ve*Time);
Ax0=Ax;
Ay0=Ay;
end
GJ=[Result(:,1),15000-Result(:,2)];
%%
clc;close all;
Ve=2940;%比冲
g=1.633; %月球重力加速度
h=600;%该阶段的下落距离
t=0; %初始时间
T=0.1; %时间步长
M_temp=[];
V_temp=[];
shijian=[];
X_temp=[];
lisan=1500:100:7500;
%lisan=5000:5100;
for i = lisan
F=i; %推力
%主减速阶段的末状态量作为快速调整阶段的初状态量
theta=55.6708*pi/180;%初速度与水平面的夹角
Vx0=32.23327;%水平初速度
Vy0=47.2005; %竖直初速度
m0=1325.255;%初始质量
Ay0=g-F*sin(theta)/(m0-F/Ve*t);%竖直初加速度
Ax0=-F*cos(theta)/(m0-F/Ve*t);%水平初加速度
count=0; %计数器
X_res=Vx0*t+0.5*Ax0*t^2;
Y_res=Vy0*t+0.5*Ax0*t^2;
Result=[];%% 迭代求 分解速度和分解位移
while (Y_res<h )
count=count+1;
Vx=Vx0+Ax0*T;
Vy=Vy0+Ay0*T;
Vx0=Vx;
Vy0=Vy;
X=Vx0*T+0.5*Ax0*T^2;
Y=Vy0*T+0.5*Ay0*T^2;
X_res=X_res+X;
Y_res=Y_res+Y;Time=count*T;
SIN=Vy/sqrt(Vy^2+Vx^2);
COS=Vx/sqrt(Vy^2+Vx^2);
Ay=g-F*SIN/(m0-F/Ve*Time);
Ax=-F*COS/(m0-F/Ve*Time);
Ax0=Ax;
Ay0=Ay;
end
M=m0-F/Ve*Time;%该阶段的末质量。
X_res; %水平位移Time=count*T; %运动时间V_res=sqrt(Vx^2+Vy^2);%合速度jiaodu=atan(Vy/Vx)*180/pi; %末速度角度%Vx %水平速度
M_temp=[M_temp;F/Ve*Time];
V_temp=[V_temp;V_res,Vx];
shijian=[shijian;Time];%记录运行时间
X_temp=[X_temp;X_res];
end
Answer=[lisan',M_temp,V_temp,shijian,X_temp];%结果总结在这里
M_temp=[];
V_temp=[];
shijian=[];
F=5085; %推力
%主减速阶段的末状态量作为快速调整阶段的初状态量
theta=55.6708*pi/180;%初速度与水平面的夹角
Vx0=32.23327;%水平初速度
Vy0=47.2005; %竖直初速度
m0=1325.255; %初始质量
Ay0=g-F*sin(theta)/(m0-F/Ve*t);%竖直初加速度
Ax0=-F*cos(theta)/(m0-F/Ve*t); %水平初加速度
count=0; %计数器
X_res=Vx0*t+0.5*Ax0*t^2;
Y_res=Vy0*t+0.5*Ax0*t^2;
Result=[];
G=[];
%% 迭代求 分解速度和分解位移
while (Y_res<h )
count=count+1;
Vx=Vx0+Ax0*T;
Vy=Vy0+Ay0*T;
Vx0=Vx;
Vy0=Vy;
X=Vx0*T+0.5*Ax0*T^2;
Y=Vy0*T+0.5*Ay0*T^2;
X_res=X_res+X;
Y_res=Y_res+Y;Time=count*T;
G=[G;X_res,Y_res,V_res,Time];
SIN=Vy/sqrt(Vy^2+Vx^2);
COS=Vx/sqrt(Vy^2+Vx^2);
Ay=g-F*SIN/(m0-F/Ve*Time);
Ax=-F*COS/(m0-F/Ve*Time);
Ax0=Ax;
Ay0=Ay;
end
GJ2=[G(:,1)+GJ(end,1),3000-G(:,2)];
flag1=[377092.173014911,2999.65005894814];%主减速段结束点
flag2=[377381.089941517,2399.99226916603];%快速调整段结束点
GJ3=[GJ;GJ2];
%plot(GJ3(:,1),GJ3(:,2));
hold on;
%plot(flag1(1,1),flag1(1,2),'o','MarkerSize',20);
%plot(flag2(1,1),flag2(1,2),'o','MarkerSize',20);
ttm=2399:-1:100;
CUBI=[zeros(2300,1)+flag2(1,1),ttm'];
GJ4=[GJ3;CUBI];
%plot(GJ4(:,1),GJ4(:,2));
flag3=[377381.089941517,100];%粗避障段结束点
%%
h=70;
aaa=1.98;
T=sqrt(2*h/g);
t1=0.5*T;
t2=t1;
x3=[];
y3=[];
x33=[];for i=0:0.01:t1
x3=[x3;0.5*1.98*i^2];
end
temp5=x3(end,1);
vvv=aaa*t1;for i=0:0.01:t2
x33=[x33;temp5+vvv*i-0.5*aaa*i^2];
end
X3=[x3;x33];
for i=0:0.01:T
y3=[y3;100-0.5*g*i^2];
end
x5=[GJ4(end,1)+X3];
GJ5=[GJ4;x5,y3];
%plot(GJ5(:,1),GJ5(:,2));
ttn=29:-1:0;
TT5=[zeros(30,1)+GJ5(end,1),ttn'];
GJ6=[GJ5;TT5];
plot(GJ6(:,1),GJ6(:,2),'LineWidth',2);
title('后4个过程的轨迹图','FontSize',14);
xlabel('水平位移/米');
ylabel('高度/米');
%axis([377092.173014911,377500,0 3000])
%% 主减速阶段的数值近似求解（用于第一问）
clc;clear;close all;
g=1.633; %月球重力加速度
m0=2.4*10^3;%卫星初始质量
Ve=2940; %比冲
theta=9.654*pi/180;%初速度与水平方向的夹角
F=7500; %推力
V0=1692.464; %近日点初速度
t=0; %初始时间
T=0.1; %时间步长
Vx0=V0*cos(-theta); %水平初速度
Vy0=V0*sin(-theta); %竖直初速度
Ay0=g-F*sin(-theta)/(m0-F/Ve*t);%竖直初加速度
Ax0=-F*cos(-theta)/(m0-F/Ve*t);%水平初加速度
count=0;
X_res=Vx0*t+0.5*Ax0*t^2;
Y_res=Vy0*t+0.5*Ax0*t^2;
Result=[];
h=12000;%主减速阶段的下落距离
%% 迭代求 分解速度和分解位移
while (Y_res<h )
count=count+1;
Vx=Vx0+Ax0*T;
Vy=Vy0+Ay0*T;
V_res=sqrt(Vx^2+Vy^2);
Vx0=Vx;
Vy0=Vy;
X=Vx0*T+0.5*Ax0*T^2;
Y=Vy0*T+0.5*Ay0*T^2;
X_res=X_res+X;
Y_res=Y_res+Y;Time=count*T;Result=[Result;X_res,Y_res,V_res,Time];
SIN=Vy/sqrt(Vy^2+Vx^2);
COS=Vx/sqrt(Vy^2+Vx^2);
Ay=g-F*SIN/(m0-F/Ve*Time);
Ax=-F*COS/(m0-F/Ve*Time);
Ax0=Ax;
Ay0=Ay;
end
M=m0-F/Ve*Time;%该阶段的末质量。
X_res %水平位移Time=count*T %运动时间V_res=sqrt(Vx^2+Vy^2) %合速度jiaodu=atan(Vy/Vx)*180/pi %末速度角度。F/Ve*Time%燃料消耗量
%% 作图
plot(Result(:,1),15000-Result(:,2),'k','LineWidth',2);
title('主减速段嫦娥三号运动轨迹图','FontSize',15);
xlabel('X轴/(m)');
ylabel('Y轴/(m)');
%% 经纬度的近似计算
R=1737.013; %平均半径
weidu=44.12-X_res/((2*pi*R)/4/90)/1000 %近月点纬度
pi*R/180; % 30.3165
%% 主减速阶段的灵敏度分析1。
clc;clear;close all;
g=1.633; %月球重力加速度
m0=2.4*10^3;%卫星初始质量
Ve=2940;%比冲
theta=9.654*pi/180;%初速度与水平方向的夹角
F=7500;%推力
RECORD=[];
lisan=11950:1:12050;
for h=lisan
%主减速阶段的下落距离h
V0=1692.464; %近日点初速度
t=0; %初始时间
T=0.1; %时间步长
Vx0=V0*cos(-theta); %水平初速度
Vy0=V0*sin(-theta); %竖直初速度
Ay0=g-F*sin(-theta)/(m0-F/Ve*t);%竖直初加速度
Ax0=-F*cos(-theta)/(m0-F/Ve*t);%水平初加速度
count=0;
X_res=Vx0*t+0.5*Ax0*t^2;
Y_res=Vy0*t+0.5*Ax0*t^2;
Result=[];%% 迭代求 分解速度和分解位移
while (Y_res<h )
count=count+1;
Vx=Vx0+Ax0*T;
Vy=Vy0+Ay0*T;
V_res=sqrt(Vx^2+Vy^2);
Vx0=Vx;
Vy0=Vy;
X=Vx0*T+0.5*Ax0*T^2;
Y=Vy0*T+0.5*Ay0*T^2;
X_res=X_res+X;
Y_res=Y_res+Y;Time=count*T;Result=[Result;X_res,Y_res,V_res,Time];
SIN=Vy/sqrt(Vy^2+Vx^2);
COS=Vx/sqrt(Vy^2+Vx^2);
Ay=g-F*SIN/(m0-F/Ve*Time);
Ax=-F*COS/(m0-F/Ve*Time);
Ax0=Ax;
Ay0=Ay;
end
M=m0-F/Ve*Time;%该阶段的末质量。
%X_res %水平位移Time=count*T; %运动时间V_res=sqrt(Vx^2+Vy^2) ;%合速度jiaodu=atan(Vy/Vx)*180/pi;%末速度角度。consume=F/Ve*Time;%燃料消耗量RECORD=[RECORD;h,consume,M,Time,V_res,X_res,jiaodu];
end
flag=320;
str={'主减速燃料消耗关于h变化图','主减速阶段末质量关于h变化图','主减速阶段总时间关于h变化图','主减速阶段末速度关于h变化图','主减速阶段水平位移关于h变化图','末速度alpha关于h变化图'};
for i=1:6flag=flag+1;subplot(flag);plot(lisan,RECORD(:,i+1),'*');title(str(i));
end
%% 各阶段时间相对误差
Fact=[487 16 125 22 19]; %实际值
Value=[421.3 24.4 62.1 9.26 3.44];%模型值
cha=abs(Fact-Value);
XD=cha./Fact;
