2014高教社杯数学建模A题解法(共27页)

1、精选优质文档-倾情为你奉上摘要本文针对嫦娥三号软着陆轨道设计与控制策略的实际问题，以理论力学（万有引力、开普勒定律、万能守恒定律等）和卫星力学知识为理论基础，结合微分方程和微元法，借助MATLAB软件解决了题目所要求解的问题。针对问题（1），在合理的假设基础上，利用物理理论知识、解析几何知识和微元法，分析并求解出近月点和远月点的位置，即139.1097 。再运用能量守恒定律和相关数据，计算出速度v1（近月点的速度）=1750.78m/s,v2（远月点的速度）=1669.77m/s，最后利用曲线的切线方程，代入点（近月点与远月点）的坐标求值，计算出方向余弦即为相应的速度方向。针对问题（2）关键词

2、：模糊评判，聚类分析，流体交通量，排队论，多元非线性回归一、问题重述嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m（见附件1）。嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现

3、软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段（见附2），要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。根据上述的基本要求，请你们建立数学模型解决下面的问题：（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。（3）对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。二、问题分析2.1问题（1）的分析首先根据问题的假设、题目中所提供的数据及图片分析

4、，可以知道嫦娥三号绕月球的轨道是由圆形轨道变为椭圆形轨道，借助开普勒定律、能量守恒定律求解出近月点的速度。为了确定近月点和元月点的精确位置及相应的速度方向，我们建立以赤道（月球的赤道）平面为xoy平面、月心为原点、月心与零度经线和零度纬线交线的交点的连线为坐标轴的坐标系和赤道（月球的赤道）平面为xoy平面，为极轴（月球的极轴）为z轴建立空间直角坐标系，x轴与极坐标系的轴相重合。首先根据着陆点的经度、纬度及月球的半径求解出着陆点和近月点（带参数a)的空间直角坐标。其次利用两点间的距离公式，并借助MATLAB软件求解出近月点与着陆点最短距离。从而计算出a（近月点的经度）=。最后利用卫星的轨迹是以月

5、心为其中一个焦点，以近月点与远月点的距离为长轴的椭圆，从而求解出卫星的轨迹方程，再运用隐函数求导的应用的知识，求解出在近月点和远月点的方向导数，进而求解近月点和远月点方向余即为近月点和远月点的速度的方向。2.2问题（2）的分析首先在根据题意，将嫦娥三号软着陆问题，分为6个阶段依次为主减速、快速调整、粗避障、精避障、缓慢下降、自由下降，我们先将6个阶段分为4个阶段，依次为第一阶段（主减速和快速调整）、第二阶段（粗避障）第三阶段(精避障),第四阶段（缓慢下降和自由下降）。其次在第一阶段粗避障阶段，嫦娥三号悬停在月球表面约2400米上方，对星下月表进行二维和三维成像，利用遗传算法的思想，从图像中先随

6、机选取部分点，能直接从三维图像中得知该点的海拔高度，再分别扫描这些点附近的地貌，找出一些地势平坦的区域，我们用区域内所有点与中心点海拔的均方差作为地势判断依据之一，保留这些坐标， 并进行重新组合，并改变某些坐标以便能获得其他新区域的坐标，再次搜索地势平坦的区域，重复进行多次搜索，直到没有出现崎岖地势的时候， 我们将此时地势最平坦的地方作为全局最优降落地点三、模型假设1、不考虑空间飞行器上各点因燃料消耗而产生的位移；2、在对卫星和空间飞行器进行轨道估计时，认为作用于其上的所有外力都通过其质心；3、卫星和空间飞行器的运动是在真空中进行的；4、卫星只受重力影响，空间飞行器除自身推力外只受重力影响；5

7、、卫星的观测图片及数据精准；6、四、变量与符号说明C0 一条车道的基本通行能力连续车流的车头间距n 条车道的基本通行能力排队长度车流量横断面通行能力系数车流量持续时间 L C y x1 x2 x3五、模型建立与求解5.1 问题(1)的分析、模型建立与求解5.1.1建模准备（1）开普勒定律开普勒第一定律开普勒第一定律开普勒第一定律，也称椭圆定律：每一个行星都沿各自的椭圆轨道环绕太阳，而太阳则处在椭圆的一个焦点中。开普勒第二定律开普勒定律开普勒第二定律，也称面积定律：在相等时间内，太阳和运动着的行星的连线所扫过的面积都是相等的。 这一定律实际揭示了行星绕太阳公转的角动量守恒。用公式表示为开普勒定律

8、开普勒第三定律开普勒定律开普勒第三定律，也称调和定律：各个行星绕太阳公转周期的平方和它们的椭圆轨道的半长轴的立方成正比。由这一定律不难导出：行星与太阳之间的引力与半径的平方成反比。这是牛顿的万有引力定a3律的一个重要基础。用公式表示为2=K开普勒定律 T这里，是行星公转轨道半长轴，是行星公转周期，是常数 。（2）万有引力万有引力：任意两个质点有通过连心线方向上的力相互吸引。该引力大小与它们质量的乘积成正比与它们距离的平方成反比，与两物体的化学组成和其间介质种类无关。即：M1M2， r2-11 其中M1，M2为两物体的质量，G=6.6710Nm.2kg.2（牛顿每平方米二次方千F=G克）5.1.

9、2 模型的建立根据以上的分析，建立以月球赤道平面为xOy平面,月心为原点O、Ox为月心与零度经线和零度纬线交线的交点的连线，Oz为极轴（月球的极轴），Oy与Ox和Oz满足右手标架，建立空间直角坐标系(如图5-1所示）。图5-1 卫星绕月轨迹及软着陆轨迹由于着陆点在球面上且近月点与远月点是由月球的经度、纬度及高度唯一确定，在此为了便于计算 将极坐标转化为空间直角坐标，并代数题中相关数据，反解出经度a。极坐标转化为空间直角坐标x=rsinjcosq即：y=rsinjsinqz=rcosj （5.1.1）x=rsin(90-b）cos(-a)y=rsin(90-b）cos(-a) （5.1.2） z

10、=rcos(90-b）距离公式：d= （5.1.3） 其中：b为纬度；a为经度；r为嫦娥三号距月心的距离；d为嫦娥三号距着陆点的距离；根据能量守恒、开普勒第二定律（面积定律），建立以下模型即： r1v1=r2v2 （5.1.4） 1122mv1+mgh=mv2+mgH22则近月点的速度，近月点的速度：v1= （5.1.5） v=2其中：m为卫星的质量，h1为海拔高度，h近月点距月球表面的距离； r1=h+r0+h1，r2=H+r0+h1，r0月球半径，H远月点距月球表面的距离， g月球重力加速度，v1 近月点的速度，v2 近月点的速度。5.1.3模型的求解5.1.3.1 近月点与远月点的位置根

11、据题目所给数据以上分析，可知：b=0,h=15000m,r0=m,h1=-2641m将以上数据代入（5.1.1）式可得，着陆点及近月点的空间直角坐标分别为： x0=r0sin(90-b)cos(-a)=r0sin(90-19.51)cos(-44.12)y0=r0sin(90-b)sin(-a)=r0sin(90-19.51)sin(-44.12) （5.1.6） z=rcos(90-b)=r0cos(90-19.51)00x=rsin(90-b）cos(-a)=(r0+h)cosay=rsin(90-b）sin(-a)=-(r0+h)sina z=rcos(90-b）=0 （5.1.7）再将

12、（5.1.6）式和（5.1.7）式代入（5.1.3）式可得关于a与d（近月点和着陆点距离）的函数，？利用Mathematica 5.0编程求解可得：a=-139.1075.1.3.2近月点与远月点的速度大小及方向近月点与远月点的速度方向，即为相应速度在x轴与y轴方向上的投影（如图5-2所示）图5-2 近月点与远月点的速度方向示意图由图易知：5.2 模型二的建立5.2.1模型准备5.2.1.1系统模型1、着陆器的动力下降段一般从15km左右的轨道高度开始，下降到月球表面的时间比较短，在几百秒范围内，所以可以不考虑月球引力摄动。月球自转速度比较小，也可忽略。因此，可以利用二体模型描述系统的运动。建

13、立图5-2所示的着陆坐标系，并假设着陆轨道在纵向平面内，令月心为坐标原点，Oy指向动力下降段的开始制动点，Ox 指向着陆器的开始运动方向。则着陆器的质心动力学方程可描述如下:r=vv=(F/m)siny-m/r2+rw2q=w w=-(F/m)cosy+2vw/r m=-F/ISP式中:r,q,w和m分别为着陆器的月心距、极角、角速度和质量；v为着陆器沿r 方向上的速度；F为制动发动机的推力(固定的常值或0）；ISP为其比m为月球引力常数；y为发动机推力与当地水平线的夹角即推力方向角。冲；图5-3 月球软着陆坐标系动力下降的初始条件由霍曼变轨后的椭圆轨道近月点确定，终端条件为着陆器在月面实现软

14、着陆。令初始时刻t0=0，终端时刻tf不定，则相应的初始条件为r0终端约束为rf=rL,vf=0,wf=0 =rL+h0,v0=0,w0=wo 式中:rL为月球半径；h0为初始轨道高度；wo为轨道角速度。月球软着陆的最优轨道设计就是要在满足上述初始条件和终端约束的前提下，调整推力大小和方向9使得着陆器实现燃料最优软着陆，即要求以下性能指标达最大。J=mdt 0tf5.2.1.2模型归一化在轨道优化过程中，由于各状态变量的量级相差较大，寻优过程中可能会导致有效位数的丢失。通过归一化处理可以克服这一缺点9，提高。计算精度。令rref=r0,mtef=m0，则=r/rref,=v/vref,vref

15、=ISp=I72=F/Fref,Fref=mrefvref/rref,=m/mref,=t/tref,=rref/vref,=q。那么，着陆器的动力学方程可改为：=v22=(F/m)siny-m/+ =w=-(F/)cosy+2/=-F/ISP相应的初始条件和终端约束变为：=1,=0,=w 000/mf=r1/r0,vf=0,wf=0性能指标改写为：=0第4期朱建丰等:基于自适应模拟退火遗传算法的月球软着陆轨道优化道优化问题转化为多参数优化问题,再利用SQP方法求解。虽然避开了没有明确物理意义的参数猜测,但是SQP的本质仍然会使该方法遇到病态梯度、初始点敏感和局部收敛问题。曾国强6和徐敏7分别

16、用二进制和浮点数GA对着陆轨道进行了优化,避免了初值猜测,得到的结果也比较满意。但是,鉴于GA局部搜索能力较差的缺点,会使得GA的优化精度不够或优化效率不高。相对而言,国外对月球软着陆轨道的优化问题研究比较少。GA最早是由Holland教授提出的8,它是一种随机优化方法,具有不依赖问题模型、适用面广和鲁棒性强的优点,并已应用在航天器的轨道优化设计中1。然而,GA在实际应用中存在收敛速度慢和早熟等问题,不具备“爬山”的能力。模拟退火算法(SAA)最早是由Kirkpatrick等提出的,它是一种启发式随机搜索算法,具有很强的局部搜索能力和“爬山”能力,但是SAA产生的新解不及GA丰富,对全局的了解

17、甚少,寻优过程很慢。因此,可以将GA和SAA的优点结合起来,扬长避短,构成高效、鲁棒的新算法。本文将GA和SAA有机地结合,形成自适应模拟退火遗传算法(ASAGA),并将其应用到月球软着陆的最优轨道设计中。1系统模型着陆器的动力下降段一般从15 km左右的轨道高度开始,下降到月球表面的时间比较短,在几百秒范围内,所以可以不考虑月球引力摄动。月球自转速度比较小,也可忽略。因此,可以利用二体模型描述系统的运动。建立图1所示的着陆坐标系,并假设着陆轨道在纵向平面内,令月心O为坐标原点,Oy指向动力下降段的开始制动点,Ox指向着陆器的开始运动方向。则着陆器的质心动力学方程可描述如下:r= vv=(F

18、/m)sin - /r2+ r 2= = -(F /m)cos+ 2v /rm= -F /ISP(1)式中:r,和m分别为着陆器的月心距、极角、角速度和质量;v为着陆器沿r方向上的速度;F为制动发动机的推力(固定的常值或0);ISP为其比冲;为月球引力常数;为发动机推力与当地水平线的夹角即推力方向角。图1月球软着陆极坐标系Fig. 1Polar coordinate system of lunar soft landing动力下降的初始条件由霍曼变轨后的椭圆轨道近月点确定,终端条件为着陆器在月面实现软着陆。令初始时刻t0= 0,终端时刻tf不定,则相应的初始条件为r0= rL+ h0,v0=

19、0,0= o(2)终端约束为rf= rL,vf= 0,f= 0 (3)式中:rL为月球半径;h0为初始轨道高度;o为轨道角速度。月球软着陆的最优轨道设计就是要在满足上述初始条件和终端约束的前提下,调整推力大小和方向,使得着陆器实现燃料最优软着陆,即要求以下性能指标达最大。J=tf0mdt (4)2归一化在轨道优化过程中,由于各状态变量的量级相差较大,寻优过程中可能会导致有效位数的丢失。通过归一化处理可以克服这一缺点9,提高计算精度。令rref= r0,mref= m0,则r= r /rref, v=v /vref,vref= /rref, ISP= ISPrref/, F= F /Fref,F

20、ref= mrefv2ref/rref, m= m /mref, = r3ref/,t= t /tref,tref= rref/vref,=。那么,着陆器的动力学方程可改写为r= vv=( F / m)sin -1 /r2+r 2= = -( F / m)cos+ 2 v /rm= - F / ISP(5)相应的初始条件和终端约束变航空学报第28卷4. 3算法的实现将提出的ASAGA应用到月球软着陆的轨道优化中,具体的实现步骤如下:(1)设置初始参数,包括种群规模M,最大遗传代数Tmax,退火初始温度T0,温度下降系数,最小新解接受次数Nmin,最大内循环次数Cmax,轨道离散参数n。随机产生

21、初始种群。(2)计算种群中各个个体的适应度值,记录最优个体,并采用如下方法进行适应度拉伸f= exp -(fmax- f)/Ti (22)式中f为拉伸后的适应度值。采用轮盘选择策略进行个体选择,并将选择的个体随机两两配对,按照式(12)、式(13)、式(15)和式(16)进行交叉操作,再利用式(14)、式(17)和式(18)对个体的每个参数进行变异操作。(3)令新解接受次数na= 0,内循环次数nc=0。对遗传操作后的子代个体计算适应度值,选择前no个优秀个体,并分别进行如下的模拟退火操作:计算(Ti),分别对个体的每个参数按式(20)进行解的变换,并按照Meteopolis准则来判断是否接受

22、新解为新的当前解。如果接受,则na= na+ 1;反之,na不变。如果na Nmin,且nc Cmax,则nc= nc+ 1,并返回 ;否则,跳出循环,并将群体中适应度最差的no个个体替换成退火操作后的个体,进入步骤)。(4)删除子代种群中的任意一个个体,并替换成步骤(2)记录的最优个体。(5)如果当前遗传代数T Tmax,则按式(19)进行降温,T= T+ 1,并返回步骤(2);否则结束整个优化过程。可以看到,在上述的ASAGA中,进行了适应度拉伸。这样处理的优点是:在温度高时(算法初期),适应度相近的个体产生的后代概率相近,避免了早期个别好的个体的后代充斥整个种群,造成早熟;而当温度不断下

23、降后,拉伸作用加强,使适应度相近的个体适应度差异放大,从而使得优秀个体更加突出。同时,算法采用了最优个体保留策略,使得最优个体不会被破坏,保证了算法的收敛性。只选择前no个优秀个体进行模拟退火操作,可以大大缩短优化时间,也能够使优秀个体向最优的方向进化。5仿真实例设定月球软着陆的初始条件为r0= 1 753km,v0= 0,0= 9. 65 10 -4 rad/s,0= 0,m0=600 kg。终端约束为:rf= 1 738 km,vf= 0,f= 0。月球引力常数= 4 902. 75km3 /s2,制动发动机推力F= 3 450 N,比冲ISP= 300 9. 8 m /s。令轨道离散化参

24、数n= 9,即待优化参数为10个推力方向角和1个终端时刻tf,共11个参数。由于寻优边界L和R决定了搜索空间的大小,而搜索空间小,则算法收敛速度快,反之则慢。因此需要估计L和R的值,加快收敛速度。根据齐奥尔科夫斯基公式和软着陆初始条件,终端时刻可由下式估计tf=1 -exp(Vf-V0)/ISP(ISPm0/F)(23)式中:Vf和V0分别为着陆器的终端速度大小和初始速度大小。同时由经验可知,推力方向角是逐渐增大的,所以寻优边界可以确定为L= 0000000000500(24)R= 1 2 3 4 5 6 7 8 9 10 9 700(25)ASAGA的其他参数设定为:M= 20,Tmax=3

25、00,T0= 10 000,= 0. 97,Nmin= 10,Cmax= 10,no= 2,(Ti)= 0. 5 1 -(T -1)/Tmax,1=52= 5 max(1,10(2T/Tmax -1),Rmax= 10。由于GA与SAA都是随机优化算法,所以ASAGA也具有随机性。因此为了验证该算法的有效性,文中对以上的初值条件,进行10次数值仿真,结果全部收敛,说明该算法有很好的收敛性。10次结果的均值为:终端速度大小Vf= 0. 11 m/s,终端月心距rf= 1 737 999. 92 m,燃料消耗277. 677 5kg。其中最好的结果如图3图6所示。图3推力方向角的时间历程Fig.

26、3Time histories of direction angle of thrust朱建丰等:基于自适应模拟退火遗传算法的月球软着陆轨道优化式中:A,B和A,B分别为父代和子代的个体;,为(0,r)区间上的均匀分布随机数,r 1为交叉系数,其大小可以控制交叉操作的变化范围;L,R分别为寻优参数的左右边界。如果交叉后的子代超出寻优边界,则重复进行交叉操作,假如重复操作达到最大设定次数Rmax,子代仍然不能满足边界约束,则利用式(13)进行子代边界修正。变异操作11采用如下形式C=C+ k(R -C)U(0,1)= 0C -k(C -L)U(0,1)= 1(14)式中:C和C分别为父代和子代的

27、个体;为(0,1)区间上的均匀分布随机数;k(0,1为变异系数;U(0,1)为随机产生的整数0或1。GA的交叉概率Pc与变异概率Pm对其性能影响很大。Pc越大,新个体产生的速度越快。然而,Pc过大会使得具有高适应度的个体结构很快被破坏;Pc过小会使搜索过程缓慢,以致停滞不前。如果变异概率Pm过小,就不易产生新的个体结构;Pm过大,GA就变成了纯粹的随机搜索算法。因此,必须采用自适应的方法,让交叉概率和变异概率随着适应度的变化而改变。同时,注意到交叉系数r与变异系数k的大小也会影响GA的性能,所以也采用自适应策略。本文中的Pc,Pm,r和k按照如下公式进行自适应调整:Pc=Pc1-Pc2(fc-

28、 favg)/(fmax-favg) fc favgPc1fc favg(15)r=r1+ r2(fmax- fc)/(fmax- favg) fc favgr1+ r2fc favg(16)Pm=Pm1-Pm2(fm-favg)/(fmax-favg) fm favgPm1fm favg(17)k=k1+ k2(fmax- fm)/(fmax- favg) fm favgk1+ k2fm favg(18)式中:fmax为群体中最大的适应度值; favg为群体的平均适应度值;fc为要交叉的两个个体中较大的适应度值; fm为要变异个体的适应度值。Pc1= 0. 9,Pc2= 0. 2,r1= 0

29、. 5 -0. 25T/Tmax(T为当前遗传代数,Tmax为最大遗传代数),r2= 0. 5,Pm1= 0. 01,Pm2= 0. 005,k1= 0. 5 -0. 4T/Tmax,k2= 0. 5。采取如上的自适应策略后,使得适应度高于群体平均适应度的个体,对应较小的r与k和较低的Pc与Pm。前者使该个体在原值附近小范围内进行交叉和变异,加快GA的收敛速度;后者使该个体得以保护进入下一代;而低于平均适应度的个体,对应较大的r与k和较高的Pc和Pm。前者使该个体进行交叉和变异的范围变大,加强GA的全局搜索能力;后者使该个体很快被淘汰掉。同时,随着优化的进行,r和k的值逐渐减小,GA的局部搜索

30、能力也逐渐加强。自适应策略使得GA在保持群体多样性的同时,保证了算法的收敛性,也在一定程度上提高了局部搜索能力。4. 2模拟退火算法SAA是基于金属退火的机理而建立起的一种随机算法,它能够通过控制温度的变化过程来实现大范围的粗略搜索与局部的精细搜索。本文采用指数降温策略对温度的变化进行控制,即Ti= T0()T-1 (19)式中:Ti为当前控制温度;T0为初始设定温度;为温度下降系数。在SAA进行时,解的变换即新解的产生,是发生在当前解的邻域内,采用如下公式进行解的变换:Y=Y+(R -Y)(Ti)U(0,1)= 0Y -(Y -L)(Ti)U(0,1)= 1(20)式中:Y和Y分别为当前解和

31、新解,为种群的个体;(Ti)为随Ti减小而减小的扰动值,当Ti为T0时,(T0) 1,保证了最大扰动值不会使新解Y越界,当Ti0时,(Ti)0,满足算法收敛的条件;为区间(0,1)上的均匀分布随机数。状态由Y变为Y的接受概率可由下面的Meteopolis准则来确定p=1 f(Y) f(Y)exp -(f(Y)-f(Y)/Ti f(Y) f(Y)(21)使用上述准则的优点是:当新解更优时,完全接受新解为新的当前解;而当新解为恶化解时,以概率p接受恶化解为新的当前解。这便使得SAA能够避免陷入局部最优。随着优化的进行,SAA的局部搜索能力也逐渐增强,确保算法有足够的搜索精度。r0= 1, v0=

32、0, 0= or30/ (6)rf= rL/r0,vf= 0,f= 0 (7)性能指标改写为J=tf0 m dt (8)3参数化方法因为月球软着陆的轨道优化搜索空间是一个泛函空间,不能直接用优化算法处理,必须将控制变量参数化。参数化方法是轨道优化与优化算法的桥梁,其精度直接影响到最优轨道的好坏程度。目前,主要的参数化方法有3种10:直接离散法、多段参数插值法和函数逼近法。直接离散法的精度最高,但计算最慢。多段参数插值法的计算最快,但精度最差。所以这两种方法都不宜采用。函数逼近法的精度适中,计算速度比较快,但是多项式系数没有明确的物理意义,很难进行初值猜测。因此,本文将函数逼近法进行改进。首先将

33、轨道离散化成许多小段,在各小段的节点设定待优化的参数,然后利用这些参数进行多项式拟合,从而得到整个轨道的控制曲线。将月球软着陆轨道离散化,分割成n个小段,每段的节点设定一个推力方向角,如图2所示。那么,可以令n+ 1个节点的推力方向角和终端时刻tf作为待优化的参数。每个节点的时刻可以由下式得到ti= t0+ i(tf-t0)/n(i= 0,1, ,n)(9)这样,就使得每个节点的推力方向角都有一个对应的节点时刻。如果假设月球软着陆的推力方向角可以表示成一个多项式,即(t)= 0+1t+2t2+3t3 (10)那么,可以将节点的推力方向角与对应的节点时刻对上述多项式进行拟合,求得多项式的系数i(

34、i= 0,1,2,3),进而得到整个着陆轨道的推力方向角。图2轨道离散化Fig. 2Trajectory discretization可以看到,采取如上改进的函数逼近法进行参数化,所选定的待优化参数都具有明确的物理意义,从而避开对没有物理含义的待优化参数的初始猜测。同时,该参数化方法的本质是函数逼近,因此它与函数逼近法的精度相当,计算速度也比较快。4自适应模拟退火遗传算法GA的局部搜索能力较差,容易陷入局部最优解,但把握搜索过程总体的能力较强;而SAA具有较强的局部搜索能力,并能使搜索过程避免陷入局部最优解,但SAA却对整个搜索空间的状况了解不多,不便于使搜索过程进入最有希望的搜索区域,从而使

35、得SAA的运算效率不高。ASAGA就是将SAA嵌入到GA的循环体中,互相取长补短,形成一种新的优化算法。4. 1自适应遗传算法GA是模拟生物在自然环境中的遗传和进化过程而形成的一种自适应全局优化随机搜索算法。由于它对搜索空间不作任何假设,它既不要求搜索空间是光滑的,也不要求它是处处可微的,因而它能解决很大一类问题,在飞行器轨道优化中逐渐开始得到应用,但是要将GA成功的应用到轨道优化中,还必须处理好如何将终端约束反映到GA的适应度函数中。考虑到轨道优化的精度,采用动态罚函数法10对终端约束进行处理,则适应度函数可表示成如下形式f(X)= -1(v(tf)-vf)2+(tf)r(tf)-frf)2

36、 1 /2 -2| r(tf)-rf|+ m0+ J(11)式中:X为种群的个体;v(tf),(tf)和r(tf)为终端时刻的状态值;1,2为惩罚因子,随着优化的进行它们逐渐变大,使得对约束的限制也逐渐加强。在遗传操作中,交叉操作是产生新个体的主要方法,它决定着GA的全局搜索能力;而变异操作只是产生新个体的辅助方法,它决定了GA的局部搜索能力。交叉与变异相互配合,共同完成对搜索空间的全局搜索和局部搜索。本文的交叉操作11采用如下方式:A=(1 -)A+BB=(1 -)B+A(12)如果A(B) R那么A(B)= R(13)第4期朱建丰等:基于自适应模拟退火遗传算法的月球软着陆轨道优化图4径向速

37、度的时间历程Fig. 4Time histories of radial velocity图5月心距的时间历程Fig. 5Time histories of distance from lunar center to lander图6角速度的时间历程Fig. 6Time history of angle velocity从结果图可以看出:虽然由ASAGA得到的最佳推力方向角与由Pontryagin极大值原理得到的最优推力方向角存在较大的差别,但是这并不影响最后的着陆效果。着陆器的终端状态为:Vf= 0. 06 m/s,rf= 1 738 000. 02 m。燃料消耗为:277. 622 4 k

38、g,比由Pontryagin极大值原理得到的最优燃料消耗(277. 576 5 kg)仅多0. 045 9kg。以上数据说明了该算法的优化精度很高。下面给出群体中最大适应度值随时间的变化曲线,如图7所示。图7说明:在优化进程中,多次出现局部“早熟”现象,但都被ASAGA的搜索能力所克服,成功地跳出了局部最优的诱惑。图7最大适应度的时间历程Fig. 7Time history of maximum fitness8.参考文献1陈刚,万自明,徐敏,等.遗传算法在航天器轨道优化中的应用J.弹道学报, 2006,18(1):1-5.Chen G, Wan Z M, Xu M, et al. Overv

39、iew of spacecrafttrajectory optimization using genetic algorithm J. Jour-nal of Ballistics, 2006,18(1):1-5.(in Chinese)2王大轶,李铁寿,马兴瑞.月球最优软着陆两点边值问题的数值解法J.航天控制, 2000(3):44-49.Wang D Y, Li T S, Ma X R. Numerical solution of TPB-VP in optimal lunar soft landing J. Aerospace Control,2000(3):44-49.(in Chin

40、ese)3王劼,崔乃刚,刘暾,等.定常推力登月飞行器最优软着陆轨道研究J.高技术通讯, 2003,13(4):39-42.Wang J, Cui N G, Liu T, et al. Study on soft-landing trajec-tories of constant-thrust-amplitude lunar probe J. HighTechnology Letters, 2003,13(4):39-42.(in Chinese)4王劼,李俊峰,崔乃刚,等.登月飞行器软着陆轨道的遗传算法优化 J.清华大学学报(自然科学版), 2003,43(8):1056-1059.Wang

41、J, Li J F,Cui N G, et al. Genetic algorithm optimi-zation of lunar probe soft-landing trajectories J. Journalof Tsinghua University(Sci and Tech), 2003,43(8):1056-1059.(in Chinese)9.附件附件1/求 m 阶B 样条：function y=N(m,x)s=factorial(m-1);sum=0;for k=0:mif(x-k)0 sum=sum+(-1)k*(nchoosek(n,k)*(x-k)(m-1); end

42、 y=s(-1)*sum/画 m 阶 B 样条：function null=l(m)x=-1:0.01:m+1;a=size(x);new=;for k=1:size(x,2) new(k)=psi(m,x(k); end plot(x,new)/对 m 阶 B 样条求 n 阶导数：function z=derive(m,n,x)sum=0; for k=0:m sum=sum+(-1)k*nchoosek(n,k)*N(m-n,x-k); end z=sum;/求 m 阶 B 样条小波:function g=psi(m,x)if m=1 if x0&x0.5&x1 g=-1;else g=0

43、; endelse s=2(-m+1);sum=0;for j=0:(2*m-2) sum=sum+(-1)j*N(2*m,j+1)*derive(2*m,m,2*x-j); end g=s*sum; end/画 m 阶样条小波的图像：function null=h(m)x=-1:0.01:2*m; a=size(x);new=;for k=1:size(x,2) new(k)=psi(m,x(k); end plot(x,new)附件2clc;clear;%绘照片图Fujian_3 = imread(附件3 距2400m处的数字高程图.tif);Fujian_4 = imread(附件4 距

44、月面100m处的数字高程图.tif);x3=1:2300;y3=1:2300;Fujian_31=double(Fujian_3);Fujian_41=double(Fujian_4);xx3 yy3=meshgrid(x3,y3);x4=0.1:0.1:100;y4=0.1:0.1:100;xx4 yy4=meshgrid(x4,y4);figure(1): mesh(xx3,yy3,Fujian_31);figure(2): mesh(xx4,yy4,Fujian_41);附件3function P = mutation(P,Nb,Pm)%变异Nbb = length(Nb);for n

45、= 1:size(P,1) b2 = 0; for m = 1:Nbb if rand Pm b1 = b2 + 1; bi = b1 + mod(floor(rand*Nb(m),Nb(m); b2 = b2 + Nb(m); P(n,bi) = P(n,bi); end endendfunction xo,fo = Opt_Simu(f,x0,l,u,kmax,q,TolFun)% 模拟退火算法求函数 f(x)的最小值点， 且 l = x = u% f为待求函数，x0为初值点，l，u分别为搜索区间的上下限，kmax为最大迭代次数% q为退火因子，TolFun为函数容许误差%算法第一步根据输

46、入变量数，将某些量设为缺省值if nargin 7 TolFun = 1e-8;endif nargin 6 q = 1;endif nargin 5 kmax = 100;end%算法第二步，求解一些基本变量N = length(x0); %自变量维数x = x0;fx = feval(f,x); %函数在初始点x0处的函数值xo = x;fo = fx;%算法第三步，进行迭代计算，找出近似全局最小点for k =0:kmax Ti = (k/kmax)q; mu = 10(Ti*100); % 计算mu dx = Mu_Inv(2*rand(size(x)-1,mu).*(u - l);%步长dx x1 = x + dx; %下一个估计点 x1 = (x1 l).*l +(l = x1).*(x1 = u).*x1 +(u x1).*u; %将x1限定在区间l,u上 fx1 = feval(f,x1); df = fx1- fx; if df 0|rand exp(-Ti*df/(abs(fx) +