高教社杯数学建模A题解法Word文档下载推荐.doc
- 文档编号:6872048
- 上传时间:2023-05-07
- 格式:DOC
- 页数:27
- 大小:103KB
高教社杯数学建模A题解法Word文档下载推荐.doc
《高教社杯数学建模A题解法Word文档下载推荐.doc》由会员分享,可在线阅读,更多相关《高教社杯数学建模A题解法Word文档下载推荐.doc(27页珍藏版)》请在冰点文库上搜索。
持续时间LCyx1x2x3
五、模型建立与求解
5.1问题
(1)的分析、模型建立与求解
5.1.1建模准备
(1)开普勒定律
开普勒第一定律开普勒第一定律开普勒第一定律,也称椭圆定律:
每一个行星都沿各自的椭圆轨道环绕太阳,而太阳则处在椭圆的一个焦点中。
开普勒第二定律开普勒定律开普勒第二定律,也称面积定律:
在相等时间内,太阳和运动着的行星的连线所扫过的面积都是相等的。
这一定律实际揭示了行星绕太阳公转的角动量守恒。
用公式表示为开普勒定律开普勒第
三定律开普勒定律开普勒第三定律,也称调和定律:
各个行星绕太阳公转周期的平方和它们的椭圆轨道的半长轴的立方成正比。
由这一定律不难导出:
行星与太阳之间的引力与半径的平方成反比。
这是牛顿的万有引力定
a3律的一个重要基础。
用公式表示为2=K开普勒定律T
这里,是行星公转轨道半长轴,是行星公转周期,是常数。
(2)万有引力
万有引力:
任意两个质点有通过连心线方向上的力相互吸引。
该引力大小与它们质量的乘积成正比与它们距离的平方成反比,与两物体的化学组成和其间介质种类无关。
即:
M1M2,r2
-11其中M1,M2为两物体的质量,G=6.67´
10Nm.2kg.2(牛顿每平方米二次方千F=G
克)
5.1.2模型的建立
根据以上的分析,建立以月球赤道平面为xOy平面,月心为原点O、Ox为月心与零度经线和零度纬线交线的交点的连线,Oz为极轴(月球的极轴),Oy与Ox和Oz满足右手标架,建立空间直角坐标系(如图5-1所示)。
图5-1卫星绕月轨迹及软着陆轨迹
由于着陆点在球面上且近月点与远月点是由月球的经度、纬度及高度唯一确定,在此为了便于计算将极坐标转化为空间直角坐标,并代数题中相关数据,反解出经度a。
极坐标转化为空间直角坐标
ì
x=rsinjcosqï
í
y=rsinjsinq
ï
z=rcosjî
(5.1.1)
x'
=rsin(90-b)cos(-a)ï
'
y=rsin(90-b)cos(-a)(5.1.2)ï
z'
=rcos(90-b)î
距离公式:
d=(5.1.3)其中:
b为纬度;
a为经度;
r为嫦娥三号距月心的距离;
d为嫦娥三号距着陆点的距离;
根据能量守恒、开普勒第二定律(面积定律),建立以下模型
r1v1=r2v2ì
(5.1.4)í
1122mv1+mgh=mv2+mgHï
î
22
则近月点的速度,近月点的速度:
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近月点与远月点的位置
根据题目所给数据以上分析,可知:
b=0,h=15000m,r0=1737013m,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)ï
00
=rsin(90-b)cos(-a)=(r0+h)cosaï
y=rsin(90-b)sin(-a)=-(r0+h)sinaï
=rcos(90-b)=0î
(5.1.7)
再将(5.1.6)式和(5.1.7)式代入(5.1.3)式可得关于a与d(近月点和着陆点距离)的函数,?
利用Mathematica5.0编程求解可得:
a=-139.107
5.1.3.2近月点与远月点的速度大小及方向
近月点与远月点的速度方向,即为相应速度在x轴与y轴方向上的投影(如图5-2所示)
图5-2近月点与远月点的速度方向示意图
由图易知:
5.2模型二的建立
5.2.1模型准备
5.2.1.1系统模型
1、着陆器的动力下降段一般从15km左右的轨道高度开始,下降到月球表面的时间比较短,在几百秒范围内,所以可以不考虑月球引力摄动。
月球自转速度比较小,也可忽略。
因此,可以利用二体模型描述系统的运动。
建立图5-2所示的着陆坐标系,并假设着陆轨道在纵向平面内,令月心为坐标原点,Oy指向动力下降段的开始制动点,Ox指向着陆器的开始运动方向。
则着陆器的质心动力学方程可描述如下:
r=vï
v=(F/m)siny-m/r2+rw2ï
q=wï
w=-[(F/m)cosy+2vw]/r⑴
m=-F/ISP
式中:
r,q,w和m分别为着陆器的月心距、极角、角速度和质量;
v为着陆器沿r方向上的速度;
F为制动发动机的推力(固定的常值或0);
ISP为其比
m为月球引力常数;
y为发动机推力与当地水平线的夹角即推力方向角。
冲;
图5-3月球软着陆坐标系
动力下降的初始条件由霍曼变轨后的椭圆轨道近月点确定,终端条件为着陆器在月面实现软着陆。
令初始时刻t0=0,终端时刻tf不定,则相应的
初始条件为
r0
终端约束为
rf=rL,vf=0,wf=0⑶
=rL+h0,v0=0,w0=wo⑵
rL为月球半径;
h0为初始轨道高度;
wo为轨道角速度。
月球软着陆的最优轨道设计就是要在满足上述初始条件和终端约束的前提下,调整推力大小和方向9使得着陆器实现燃料最优软着陆,即要求以下性能指标达最大。
J=ò
mdt0tf
5.2.1.2模型归一化
在轨道优化过程中,由于各状态变量的量级相差较大,寻优过程中可能会导致有效位数的丢失。
通过归一化处理可以克服这一缺点[9],提高。
计算精度。
令rref=r0,mtef=
m0,则=r/rref,=v/vref,vref=ISp=I7
2=F/Fref,Fref=mrefvref/rref,=m/mref,==t/tref
=rref/vref,=q。
那么,着陆器的动力学方程可改为:
=vï
22=(F/m)siny-m/
+ï
ï
=wï
=-[(F/)cosy+2]/ï
=-F/ISP相应的初始条件和终端约束变为:
=1,=0,=w000/m
f=r1/r0,vf=0,wf=0
性能指标改写为:
=ò
第4期朱建丰等:
基于自适应模拟退火遗传算法的月球软着陆轨道优化
道优化问题转化为多参数优化问题,再利用SQP
方法求解。
虽然避开了没有明确物理意义的参数
猜测,但是SQP的本质仍然会使该方法遇到病态
梯度、初始点敏感和局部收敛问题。
曾国强[6]和徐
敏[7]分别用二进制和浮点数GA对着陆轨道进行
了优化,避免了初值猜测,得到的结果也比较满意。
但是,鉴于GA局部搜索能力较差的缺点,会使得
GA的优化精度不够或优化效率不高。
相对而言,
国外对月球软着陆轨道的优化问题研究比较少。
GA最早是由Holland教授提出的[8],它是
一种随机优化方法,具有不依赖问题模型、适用面
广和鲁棒性强的优点,并已应用在航天器的轨道
优化设计中[1]。
然而,GA在实际应用中存在收
敛速度慢和早熟等问题,不具备“爬山”的能力。
模拟退火算法(SAA)最早是由Kirkpatrick等提
出的,它是一种启发式随机搜索算法,具有很强的
局部搜索能力和“爬山”能力,但是SAA产生的
新解不及GA丰富,对全局的了解甚少,寻优过程
很慢。
因此,可以将GA和SAA的优点结合起
来,扬长避短,构成高效、鲁棒的新算法。
本文将GA和SAA有机地结合,形成自适应
模拟退火遗传算法(ASAGA),并将其应用到月
球软着陆的最优轨道设计中。
1 系统模型
着陆器的动力下降段一般从15km左右的轨
道高度开始,下降到月球表面的时间比较短,在几
百秒范围内,所以可以不考虑月球引力摄动。
月
球自转速度比较小,也可忽略。
因此,可以利用二
体模型描述系统的运动。
建立图1所示的着陆坐
标系,并假设着陆轨道在纵向平面内,令月心O
为坐标原点,Oy指向动力下降段的开始制动点,
Ox指向着陆器的开始运动方向。
则着陆器的质
心动力学方程可描述如下:
•r=v
•v=(F/m)sinψ-μ/r2+rω2
•θ=ω
•ω=-[(F/m)cosψ+2vω]/r
•m=-F/ISP
(1)
r,θ,ω和m分别为着陆器的月心距、极角、
角速度和质量;
v为着陆器沿r方向上的速度;
F
为制动发动机的推力(固定的常值或0);
ISP为其
比冲;
μ为月球引力常数;
ψ为发动机推力与当地
水平线的夹角即推力方向角。
图1 月球软着陆极坐标系
Fig.1 Polarcoordinatesystemoflunarsoftlanding
动力下降的初始条件由霍曼变轨后的椭圆轨
道近月点确定,终端条件为着陆器在月面实现软
着陆。
令初始时刻t0=0,终端时刻tf不定,则相
应的初始条件为
r0=rL+h0,v0=0,ω0=ωo
(2)
rf=rL,vf=0,ωf=0(3)
rL为月球半径;
h0为初始轨道高度;
ωo为
轨道角速度。
月球软着陆的最优轨道设计就是要在满足上
述初始条件和终端约束的前提下,调整推力大小
和方向,使得着陆器实现燃料最优软着陆,即要求
以下性能指标达最大。
J=∫tf0•mdt(4)
2 归一化
在轨道优化过程中,由于各状态变量的量级
相差较大,寻优过程中可能会导致有效位数的丢
失。
通过归一化处理可以克服这一缺点[9],提高
令rref=r0,mref=m0,则–r=r/rref,v=
v/vref,vref=μ/rref,ISP=ISPrref/μ,F=F/Fref,
Fref=mrefv2ref/rref,m=m/mref,ω=ωr3ref/μ,–t=t/
tref,tref=rref/vref,–θ=θ。
那么,着陆器的动力学方
程可改写为
–r=v
v=(F/m)sinψ-1/–r2+–rω2
–θ=ω
ω=-[(F/m)cosψ+2vω]/–r
m=-F/ISP
(5)
相应的初始条件和终端约束变
航 空 学 报第28卷
4.3 算法的实现
将提出的ASAGA应用到月球软着陆的轨
道优化中,具体的实现步骤如下:
(1)设置初始参数,包括种群规模M,最大遗
传代数Tmax,退火初始温度T0,温度下降系数κ,
最小新解接受次数Nmin,最大内循环次数Cmax,轨
道离散参数n。
随机产生初始种群。
(2)计算种群中各个个体的适应度值,记录
最优个体,并采用如下方法进行适应度拉伸
f′=exp[-(fmax-f)/Ti](22)
式中f′为拉伸后的适应度值。
采用轮盘选择策
略进行个体选择,并将选择的个体随机两两配对,
按照式(12)、式(13)、式(15)和式(16)进行交叉操
作,再利用式(14)、式(17)和式(18)对个体的每个
参数进行变异操作。
(3)令新解接受次数na=0,内循环次数nc=
0。
对遗传操作后的子代个体计算适应度值,选择
前no个优秀个体,并分别进行如下的模拟退火
操作:
①计算δ(Ti),分别对个体的每个参数按式
(20)进行解的变换,并按照Meteopolis准则来判
断是否接受新解为新的当前解。
如果接受,则
na=na+1;
反之,na不变。
②如果na<
Nmin,且nc<
Cmax,则nc=nc+1,
并返回①;
否则,跳出循环,并将群体中适应度最
差的no个个体替换成退火操作后的个体,进入步
骤)。
(4)删除子代种群中的任意一个个体,并替
换成步骤
(2)记录的最优个体。
(5)如果当前遗传代数T<
Tmax,则按式(19)
进行降温,T=T+1,并返回步骤
(2);
否则结束整
个优化过程。
可以看到,在上述的ASAGA中,进行了适
应度拉伸。
这样处理的优点是:
在温度高时(算法
初期),适应度相近的个体产生的后代概率相近,
避免了早期个别好的个体的后代充斥整个种群,
造成早熟;
而当温度不断下降后,拉伸作用加强,
使适应度相近的个体适应度差异放大,从而使得
优秀个体更加突出。
同时,算法采用了最优个体
保留策略,使得最优个体不会被破坏,保证了算法
的收敛性。
只选择前no个优秀个体进行模拟退
火操作,可以大大缩短优化时间,也能够使优秀个
体向最优的方向进化。
5 仿真实例
设定月球软着陆的初始条件为r0=1753
km,v0=0,ω0=9.65×
10-4rad/s,θ0=0,m0=
600kg。
终端约束为:
rf=1738km,vf=0,ωf=0。
月球引力常数μ=4902.75km3/s2,制动发动机
推力F=3×
450N,比冲ISP=300×
9.8m/s。
令轨道离散化参数n=9,即待优化参数为10
个推力方向角和1个终端时刻tf,共11个参数。
由于寻优边界L和R决定了搜索空间的大小,而
搜索空间小,则算法收敛速度快,反之则慢。
因此
需要估计L和R的值,加快收敛速度。
根据齐奥
尔科夫斯基公式和软着陆初始条件,终端时刻可
由下式估计
tf={1-exp[(Vf-V0)/ISP]}(ISPm0/F)
(23)
Vf和V0分别为着陆器的终端速度大小和
初始速度大小。
同时由经验可知,推力方向角是
逐渐增大的,所以寻优边界可以确定为
L=[0 0 0 0 0 0 0 0 0 0 500]
(24)
R=[[12345678910]×
9700]
(25)
ASAGA的其他参数设定为:
M=20,Tmax=
300,T0=10000,κ=0.97,Nmin=10,Cmax=10,
no=2,δ(Ti)=0.5×
[1-(T-1)/Tmax],σ1=
5σ2=5max(1,10(2T/Tmax-1)),Rmax=10。
由于GA
与SAA都是随机优化算法,所以ASAGA也具
有随机性。
因此为了验证该算法的有效性,文中
对以上的初值条件,进行10次数值仿真,结果全
部收敛,说明该算法有很好的收敛性。
10次结果
的均值为:
终端速度大小Vf=0.11m/s,终端月
心距rf=1737999.92m,燃料消耗277.6775
kg。
其中最好的结果如图3~图6所示。
图3 推力方向角的时间历程
Fig.3 Timehistoriesofdirectionangleofthrust
朱建丰等:
A,B和A′,B′分别为父代和子代的个体;
α,β为(0,r)区间上的均匀分布随机数,r≤1为交
叉系数,其大小可以控制交叉操作的变化范围;
L,R分别为寻优参数的左右边界。
如果交叉后
的子代超出寻优边界,则重复进行交叉操作,假如
重复操作达到最大设定次数Rmax,子代仍然不能
满足边界约束,则利用式(13)进行子代边界修正。
变异操作[11]采用如下形式
C′=C+kγ(R-C) U(0,1)=0
C-kγ(C-L) U(0,1)=1 (14)
C和C′分别为父代和子代的个体;
γ为
(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-favg)/(fmax-favg)fc≥favg
Pc1fc<
favg
(15)
r=r1+r2(fmax-fc)/(fmax-favg)fc≥favg
r1+r2fc<
(16)
Pm=Pm1-Pm2(fm-favg)/(fmax-favg)fm≥favg
Pm1fm<
(17)
k=k1+k2(fmax-fm)/(fmax-favg)fm≥favg
k1+k2fm<
(18)
fmax为群体中最大的适应度值;
favg为群体
的平均适应度值;
fc为要交叉的两个个体中较大
的适应度值;
fm为要变异个体的适应度值。
Pc1=0.9,Pc2=0.2,r1=0.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的局部搜索能力也逐渐加强。
自适应策
略使得GA在保持群体多样性的同时,保证了算
法的收敛性,也在一定程度上提高了局部搜索
能力。
4.2 模拟退火算法
SAA是基于金属退火的机理而建立起的一
种随机算法,它能够通过控制温度的变化过程来
实现大范围的粗略搜索与局部的精细搜索。
本文
采用指数降温策略对温度的变化进行控制,即
Ti=T0(κ)T-1(19)
Ti为当前控制温度;
T0为初始设定温度;
κ
为温度下降系数。
在SAA进行时,解的变换即
新解的产生,是发生在当前解的邻域内,采用如下
公式进行解的变换:
Y′=Y+(R-Y)δ(Ti)ξ U(0,1)=0
Y-(Y-L)δ(Ti)ξ U(0,1)=1
(20)
Y和Y′分别为当前解和新解,为种群的个
体;
δ(Ti)为随Ti减小而减小的扰动值,当Ti为
T0时,δ(T0)≤1,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 高教 数学 建模 题解