三号软着陆轨道设计与控制策略

三号轨道设计与控制策三号在高速飞行的情况下，要保证准确地在月球预定区域内实现，关键问题是着陆轨道与控制策略的设计。本文从题目所三个问题出发，建立数学模型，求解相关着陆轨道，确定最优控制策略。针对第一个问题，我们运用了模型，最后求出v近=1.9km/,v远=.61k/s然后，我们再根据探测器的降落点位置（19.51°W，44.12°N，考虑过程中飞船位置的变化，由一个相对比较粗糙的求解最后确定近月点的位置为19.4W4.4N)，远月点的位置(160.46E,5.6S)。针对第二个问题，我们着重研究了主阶段，快速阶段以及避障阶段的控制策略。主阶段通过架构在二维平面上探测器的运动规律，建立动力学模型，并通过直接法求解。轨道离散化之后，将泛函极值问题转化为多约束多变量非线性规划最值问题，通过数值计算微分方程的unge-ua法和针对本问题改进的遗传算法求解最优控制策略和轨迹。快速阶段由于路径、时间较短，所以未进行优化，只是进行了时间和能耗的计算。最终算出时间t≈21,能耗约为33g。避障阶段，首先用聚类算法预处理数据，排除了明显的“坑”和“山”，然后通过计算坡度，综合衡量安全度和区域大小，找出最适合区域O(1)–O(0.1)kmO(10)–O(1)m/s，如果进化样本足够多，可以达到设置的最优精度。而Runge-Kutta法则可将微分方程的数值计算解精确到2–4阶。而无论是遗传算法或是Runge-Kutta法都具有较高的敏感性，不依赖于初值，算法稳定。本题模型将理论与仿结合，在一定误差允许范围内得到最优方案，具有较强的可操作性。同时具有良好的稳定性和敏感性。但由于为了简化计算，牺牲了一定的精度，同时轨道也并不完善。：，轨道 遗传算法，避问题提月球是距离地球最近的天体，是人类深空探测的第一站。在各国纷纷制定月球探测计划同时，我国探月工程总体规划是在2020年前实现“绕、落、回”3个发展阶段，而三号需要完成的便是月球目标。三号在高速飞行的情况下，要保证准确地在月球预定区域内实现，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其过程共分为6个阶段，要求满足每个阶段在关键点所处的状态；尽量减少过程的消耗。根据上述基本要求，建立数学模型解决下面的问题：确定着陆准备轨道近月点和远月点的位置，以及三号相应速度的大小与方向确定三号的着陆轨道和在6个阶段的最优控制策略问题分本题第一问要求确定着陆准备轨道近月点和远月点的位置，以及三号相应速度大小和方向。相应的速度和方向可以根据动力学进行分析，这里主要应用了模型。而近月点和远月点位置需要在分析主阶段后确定。第二问在确定三号的着陆轨道和各阶段的最优控制策略中，我们进行各个阶段的展开分析，其中着重解决主阶段和避障阶段问题，其他阶段由于时间短或过程简单原因不做详细说明。误差分析和敏感性分析分别会在我们处理各个阶段的过程中进行说明。模型假假设过程中月球探测器在一个固定的铅垂面内运动，并且不考虑侧向运动和月球自转。符号说V(r)：探测器在月球势场中的势能 探测器在着陆准备轨道上的机械能 探测器在着陆准备轨道上的角动量rmin：探测器在着陆准备轨道上的近月点到月心的距离rmax：探测器在着陆准备轨道上的远月点到月心的距离v：探测器在着陆准备轨道上的近月点的速度v F 发的推力大 比冲 即单位质量的推进剂产生的推模型建立与求该问题为问题，以月球为原点建立极坐标系，近月点与月球的连线为极轴。如图示Figure1:系统的有效势的形式

Veff(r)=V(r)+2mr2将月球当作标准球形，则V(r)=α，其中的α=GMmr已知在探测器绕月球转动过程中能量、角动量分别守恒，得E=1m(˙2+r2˙2)+V 2J=mr2˙2=C(常量 对于方程（1)、(2)，消去时间并积分，得到的结果是：ppe

=1+ecos rp=J2 e

1+mα2方程 是极坐标中的标准圆锥曲线方程，极坐标原点位于圆锥曲线的焦点。反解得√ C E=GMm(e2− (rmin0)(rmax,π)(3),ep，再将求得的e代入公式(4)得到能量E。又因为：E=1mv2− rminrmaxvv远。v1.69km/sv1.61km/s方向如上图中所示，分别垂直于相应r.位置的求解将在下面阶段中给确定三号的着陆轨道和在 个阶段的最优控制策主阶一、模型本题中探测器的动力学模型采用二维模型[1]，即假设月球探测器在一个固定的铅垂面内的基础上运动，没有考虑侧向运动，而且得到的模型都是在忽略月球自传的基础上得到的。我们可以建立如下所示的坐标系：Figure2:坐o为坐标原点，oy指向着陆转移轨道的近月点：r∈R+为探测器到月心的距离；θoyor的夹角；ψ(t)为推力方向or垂线的夹角；F为发的常值推力大小 着陆器质心运动学方程为r˙=v˙=Fsinψ−µ+ rC˙=−C其中v∈Rr方向上的速度；ω∈RθmR+是探测器质量；µ是月球引力常数；C为制动火箭的排气速度，是一个常值。假定初始时刻t0=0，终端时刻tf为任意值。的初始条件由探测器在椭圆轨道近月点的状态r(0)= v(0)= θ0= ω(0)= m(0)=为了在到达月面时实现，有如下终端条r(tf)=rf v(tf)=vf ω(tf)=r0是月心到近月点的距离，rf是月球半径，vf最优轨道设计的目的是寻找最优控制u(t)=u∗(t)，其中u(t)=[ψ,tf]T，f∫f取最小二、模型

J=

˙dt=m0−m(tf月球探测器主降速阶段优化问题是一个泛函优化问题，一般可通过直接法，间接法或者两者混合法求解[2]。间接法是指利用变分法或者otrain极大值原理将轨道优化问题归为终端时间自由的两点边值问题，然后再对此求解。但是这种方法往往要猜测一个极值，而且初值并没有物理意义，这就使得计算容易陷入局部最优，或者难以收敛。而直接法是指直接对动力学方程和控制量进行离散化处理，通过对离散参数进行优化来得到最优控制方案和最优轨迹。[3]n本文采用的是直接法求解。首先，将轨道离散化，分割成 个等时段，n应的，每个小节点的推力方向为

0,·,

,∆t=tf−t0,ti=t0+ i=1···n,3于是，这个泛函极值优化问题转变为多约束多变量非线性规划问题，即已有上述模型中的初值条件和终值条件的约束，加上动力学方程组的条件约束，求解变量[ψ0,·,ψn,tf],使得性能指标函数J最大。为了求解这一问题，我们打算采用在处理此类问题中发展已经比较成算法Figure3:Figure4:参数初这里需要初始化的参数包括离散化参数N，种群规模popsize，最大遗传代数Generationnmax,交叉pcrossover,变异pmutation,precision,参数的左右界L，R。产生初始由于探测器初始切向速度很大，因此，我们不妨把一个随机产生的递增的角度作为初始种群，这样可以在一定程度上避免算法难以收敛的情况。我们根据给定的初始条件和动力学方程组，来求解系统的末状态参量，再根据末状态参量来计算种群的适应度。这里，我们的适应度函数通过性能指标函数，并

Figure加入与给定的终端状态相比得到的罚函数来共同刻画，可√fitness=J− v2+(ω∗r)2−Vf|−q2∗|r−rf而在求解动力学方程组时，为了提高精度又不致计算量太大，我们采用了二阶Runge-Kutta方法来求解非线性选择操这里，我们的选择操作采用的是基于赌法的非线性选择，以适应度值作为依据。交叉变这里交叉采用多点交叉和均匀交叉，并逐渐增加交叉变异N= popsize= Generationnmax=1000,F=pcrossover= pmutation= precision=L= R= 为了使算法快速收敛并且控制得到数据的精度，在适应度函数中加入关于rVrO(1)km，V的精度约束至的精度，但是限于计算能力以及优化范围有限，我们暂时选取比较粗糙的精度。如图5所示，可以看出，经过算法有较好的收敛性。为了克服算法在未到达指定终值条件精度范围内出现早熟现象，我们拟打算以大量实验，排除不符合终值精度的进化，最后再根据满足条件的解的适应度值大小作比较，得出最优轨迹。该算法适用于并行计算，如果有条件，那么效率将大大提升。限于条件因素，我们做了30次试验。从实验结果看出，有大量的结果不满足精 作为我们算法的分析评判依据。我们得到了如下的径向距r和切向角θ变化图

Figure6:从图6中我们可以看到，θ的角度变化比较平缓，变化速度也符合规律。但是r的轨迹有一小段波动，从方位角的设定中，我们可以看到原因，即方位角之间过于离散，没有良好的连续的变化，导致系统运动不是很稳定。但是r的精度在我们可以忍受的范围内，故也算是比较好的一个结果。在次的大量重复试验后，我们会得到一个最优的结果三、近地点位置的我们假设探测器降在同一个平面上进行的，而且其他着陆过程相对于主阶段切向位移要小很多，所以我们暂且忽略不计，于是，由主阶段得到的θ我们可以粗略地推算出探测器在着陆过程中的纬度变化。根据我们模拟出来的结果，这一纬度变化值为9.5949°。而我们通过计算出来的过程的总时间，考略到月球的自转，粗略地估计经度的变化。这里，在主降速阶段精度的变化为0.0305°。四、精度与敏感性由于遗传算法不受初值的影响，故其灵敏度较高。利用遗传算法求解的角度和时间的精度为减少计算量我们设为0.1，一般这个精度可设为1e4。但是，由于遗传算法的局部搜索能力较差，所以其精度很难保证。二阶Runge-Kutta法数值比较稳定，对初值依赖性较低，顾灵敏度较高。而利用二阶Runge-Kutta法差分求解的精度为2阶的，而其中涉及到的变量的基本单位是m，m/s，kg，s，因此，二阶的精度适合我们较为错略地估算。如果条件允许，最好采用四阶Runge-Kutta法来求解，这样精度可以精确到4阶，但四阶以上的不推荐，因为不稳定五、进一本文我们使用了遗传算法来求解这一优化问题，但事实上，遗传算法计算量大，局部搜索能力较差，且容易出现早熟现象，这对于我们更加严苛的精度来说，是不利的。对此，我们可以考虑在此基础上，对遗传算法的过程进行更进一步的优化，将算法更加有机地融合到这个过程中，更加自适应。或者也可以考虑将其他优化算法与遗传算法相结合来进行优化。比如说，可以引进模拟退火，来弥补遗传算法在局部搜索能力差的缺陷[4]，或者是引进人工神经网络算法来继续优化等等，这里不再赘述由于该阶段要迅速改变探测器的姿态使其由水平变为垂直，同时推力从最大变为最低，可假设该阶段推动力(F)和探测器俯仰角(θ)均随时间线性变化，且由于时间短暂，质量变化忽略不计。初始状态为(1,1),该阶段结束时的状态为(2,2)，设该阶段用时为t0，由线性变化{将末状态带入上方公式{

F=F1−θ=θ1−a=a(t)

对水平方向进行力学分析

b=b(t)=∫m(v2x−v1x) (F1−at)sin(θ1−0∫— t(F−at)sin(θ1−0 0 F1(cosθ2−cosθ1)−a[t0cosθ2+1(sinθ2−sin21=bF1t0(cos21

1)–1)

F1−F2

022θ t0(sin022

–sinθ

θ−θ θ−1 12 t0[F2cosθ2−F1cosθ1−(Fθ−θ θ−1 12t —

mv1x(θ1−(F1−F2)(sinθ2−sin带入以下

F2cos

1cosθ1

F1= F2= θ1≈ θ2≈10◦.由上面主部分求得m=求则相应的消耗总质量

t0≈ ∫t0F ∫t0 1(F1t0−

a≈280m/s, t0≈21s, ve=2940m/sm=33kg由该结果可知燃耗质量与总质量相比几乎可以不计。避障阶基础信息在本阶段中，先读入附件中的，用等程序处理后用聚类算法[6]（见附录三）算出可行区域，得到粗略的可行区域解以及空格分割的等高线表（in文件），并坡度的±5xy方向上的平均坡度，并将其平方和Figure突出点（山和坑）的计算和评分对于图中每个11∗11的区域，对该区域的每个点按照中心点的高度和坡论计算，与实际值相减，取绝对值累加入改点的”不平均得分“中.这样每个点被121个11∗11的覆盖率累加得分。最后控制条件（坡度小于tan25，累加得分小于区域的30%(比例不定））选出可行区域。以01点阵的形式输出可行区域，并对可行区域进行部分和计数。该输出结果可直接使用，也可以用excel或其它可视化软件转化为。(如下图9） FigureFigure9:扫描确定逐步增大正方形边长的值，并以此扫描区域，并设置适当的安全系数使得区域有合适的大小，直至找到当前安全系数下最大的可行区域。之后在所有的最大可行区域中找到距离地图中心点最近的一个以节约能量。最后得到以下结果：模型的评价与时空复杂度可改进的地方：边界的处理（这个程序里由于边界不到1%的长度或者0.01%的加lgn倍时间复杂度为代价进行圆形区域的搜索（圆形区域优于正方形区域。另外，在搜索可行区域时也可以用循环队列作为，在上一层的可行区域基础上对下一层时间复杂度分析：该程序的理论时间复杂度为O(n3)，但由于最终搜索到的可行区域不会太大，所以在正常的运行时时间复杂度一般会缩减到O(kn2)，k«n。空间复杂度分析；该程序的空间复杂度即为输入大小， 可移植性分析：在不同的平台、环境上使用时，可能需要根据目的地的大致环境，少许修改平坦度参数、凹凸得分系数、百分比等参数，以寻找的合适的降 点。Figure10:缓速下降阶段主要考虑到着陆安全性，以尽可能小的速度匀速垂直下降，一般选择速度为-2/,水平速度为0，一般只需要控制加速度稍小于月球引力加速度即可。[5]而自由落体阶段就是不加推力让探测器自由下落，这两阶段相对比较简单，我们便不再展开。模型分模型优理论与仿结合，可操作性强。本文通过物理背景建立理想化的模型，将艰涩的模型理论求解转化为可操作的数值仿真模拟计算，在一定误差允许范围内得到一个最优的结果。有良好的稳定性和敏感性。文中所涉及主轨道的遗传算法，以及避障算法都对初值无依赖性，有良好的敏感性；同时遗传算法作为发展比较成算法也具有极强的稳定性。模型缺精度较差，为了简化计算，我们这里做了扩大精度范围的牺牲，不过可以说明我们的方法适合更的精度，但需要不断优化过程及其修改罚函数系数。轨道不够完善。由于时间紧迫，我们着重考虑了预备阶段、主阶段、快速调整阶段的轨道和避障技术的实现，但并没有很好地完善剩下阶段的轨道情况，但由于其相比之下比较简单，故我们不妨暂且忽略。参考文[1],.月球轨道优化方法比较研究,工程数学学报,第29卷第期[2],.基于自适应模拟退火遗传算法的月球轨道优化,航空学报,第28卷第4期,2007.[3],.月球轨道优化方法比较研究,工程数学学报,第29卷第期[4],.基于自适应模拟退火遗传算法的月球轨道优化,航空学报,第28卷第4期,2007.[5].三号着陆器动力下降的制导导航与控制,中国科学,第44卷第期:377-[6].三号自主避障控制技术,中国科学,第44卷第6期:559-附录%月球遗传算法优%参数初始globalmnNewPopchildren1children2VarNumNGGenerationnmaxvrwthetaglobalboundsJ %离散化参 %精%参数的左右L=[000000000L=[LR=[203040506070809090R=[Rbits=ceil(log2((bounds(:,2)-bounds(:,1))'./precision));%由设定精度划分区间whileG<=Generationnmaxforj=1:mvalue(1,j)=fitness(b2f(Pop(j,:),bounds,bits));%[selectpop]=selection(Pop,bounds,bits);%非线性选[CrossOverPop]=crossover(selectpop,pcrossover,round(unidrnd(Generationnmax-pMutation=pm0+(G^4)*(pcrossover/3-ifG>Generationnmax/2&&MaxValue==0t=1:(G-plot(t,Trace(1:(G-%title('函数优化的遗传算法xlabel('进化世代数(Generationnmax)');ylabel('每一代最优适应度%holdon;%text(I+5,MaxFval,['FMAX='%str1=sprintf('进化到%d代,自变量为%s时,得本次求解的最优值%f\n对应是%figure(2);plot(t,p);%绘制变异值增大过程functionf=fitness(X)globalNGGenerationnmaxvrwthetaJ r0=15*10^3+rL;%km 待 mu=4902.75*10^9;%m3/s2Isp=2940;%m/s dm=-Vf=57;%m/s%for vv=bav+h*(baF/(bam0-badm*i*X(1,11)/N)*sin(X(1,i))-%%J=bam0-%f=bam0+J-q1*sqrt((bav-bavf)^2+(baw*bar-bawf*barf)^2)-q2*abs(bar-%f=J-q1*abs(sqrt(bav^2+(baw*bar)^2)-Vf)-q2*abs(bar-forj=1:Nforf=m0-J-q1*abs(sqrt(v^2+(w*r)^2)-Vf)-q2*abs(r-%f=q1*abs(sqrt(v^2+(w*r)^2)-Vf)+q2*abs(r-iff<0ifpara>2%&¶<ifabs(r-rf)>100*10^3ifabs(sqrt(v^2+(w*r)^2)-

if

ifabs(r-rf)>100^(para-1)*10^3ifabs(sqrt(v^2+(w*r)^2)-Vf)>100^(para-1)if

%初始化种function[initpop]=InitPop(popsize,bits)globalinitpop=zeros(popsize,len);%Thewholezeroencodingindividualfori=1:popsize-1whilemark==1forj=1:9ifXtemp(1,j)>Xtemp(1,j+1)

%数值转function[fval]=b2f(bval,bounds,bits)cs=[0cumsum(bits)];fori=1:numV%选择操%采用基于赌法的非线性选%各成员按适应值从大到小分配选择概率%P(i)=(q/1-(1-q)^n)*(1- 其中P(0)>P(1)>...>P(n),function[selectpop]=selection(pop,bounds,bits)globalmnforfit(i)=fitness(b2f(pop(i,:),bounds,bits));%以函数值为适应值做依selectprob=fit/sum(fit);%计算各相对适应度(0,1)[ywhilenewIn<=m&&ifrNums(newIn)<newfit(fitIn)

%交叉操functionglobalmiflen>2&&mod(len,2)==1%如果用来进行交叉的的条数为奇数，将其调整为偶iffori=0:2:length(y1)-2ifopts==0

%采用均匀交functionglobalchildren1(crossposition)=parent1(crossposition);%掩码为1，父1为子1提供children2(crossposition)=parent2(crossposition);%掩码为1，父2为子2提供%多点交functionglobalnVarNumfori=1:VarNum%变异操functionglobalmniffork=unidrnd(n,1,VarNum);%设置变异点数，一般设置1forifOldPop(position(i),k(j))==1fori=1:30Y=[XrvwthetaJG];fprintf(fid,'%f',Y); r0=15*10^3+rL;%km 待 %1500-dm=- %for vv=bav+h*(baF/(bam0-badm*i*X(1,11)/N)*sin(X(1,i))-%%J=bam0-%f=bam0+J-q1*sqrt((bav-bavf)^2+(baw*bar-bawf*barf)^2)-q2*abs(bar-%f=J-q1*abs(sqrt(bav^2+(baw*bar)^2)-Vf)-q2*abs(bar-forj=1:Nfor%f=J+q1*abs(sqrt(v^2+(w*r)^2)-Vf)+q2*abs(r-%f=q1*abs(sqrt(v^2+(w*r)^2)-Vf)+q2*abs(r-%用于画X=[0.0000000.0000000.6349210.7936510.0000000.00000032.75590651.02362213.46456722.677165400.391389];rL=1737.013*10^3;%m%初始条r0=15*10^3+rL;%kmw0=1693/r0; 待rf=3*10^3+rL;%mmu=4902.75*10^9;%m3/s2Isp=2940;%m/s dm=-Vf=57;%m/s%for vv=bav+h*(baF/(bam0-badm*i*X(1,11)/N)*sin(X(1,i))-%%J=bam0-%f=bam0+J-q1*sqrt((bav-bavf)^2+(baw*bar-bawf*barf)^2)-q2*abs(bar-%f=J-q1*abs(sqrt(bav^2+(baw*bar)^2)-Vf)-q2*abs(bar-Vimg=[Vimgv0];Thetaimg=[Thetaimgtheta];forj=1:Nforvv=v+h*(F/(m0-dm*i*h)*sin(X(1,j))-Vimg=[VimgRimg=[RimgThetaimg=[Thetaimgtheta];Wimg=[Wimgw];

%f=J+q1*abs(sqrt(v^2+(w*r)^2)-Vf)+q2*abs(r-%f=q1*abs(sqrt(v^2+(w*r)^2)-Vf)+q2*abs(r-附录{$inline{$R-,S-,M-r=11;//averageranger2=rshr1;sr=5;//safetyrange

w,dx,dy,s:array[-r..nx+r+1,-r..ny+r+1]oflongint;//s:sumk,q,sq:array[0..nx,0..ny]oflongint;adx,ady,b,ad:array[-r..nx+r+1,-r..ny+r+1]ofextended;//b:blockansx,ansy:array[0..nx*ny]oflongint;qu:array[1..nx*ny]ofproceduresort(l,r:longint);i,j:longint;x:=qu[(l+r)div2];whilequ[i]<xdowhilex<qu[j]doifnot(i>j)thenj:=j-untili>j;ifl<jthenifi<rthenfori:=1tonxdoforj:=1tonydo

fori:=1tonxdoforj:=1tonydoread(k[i,j]);fori:=1tonxdoforj:=1tonydobeginforii:=1tor2doadx[i,j]:=adx[i,j]+w[i+ii,j]-w[i-ady[i,j]:=ady[i,j]+w[i,j+ii]-w[i,