2014高教社杯全国大学生数学建模竞赛
承 诺 书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A
我们的报名参赛队号为（8位数字组成的编号）：
所属学校（请填写完整的全名）：
参赛队员 (打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人 (打印并签名)：
（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。）
日期： 2014 年9月 14 日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2014高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
摘要
嫦娥三号在登月的整个过程中，最关键的是其软着陆的过程，而嫦娥三号在高速飞行的情况下，保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计，因此本项目研究意义重大。
针对问题一，由于飞行器绕月球运动的过程与卫星绕地球运行的过程相似，因此我们可使用卫星绕地球运动的研究方法，地球卫星入轨后，处于地球的保守场中，在不计阻力损耗情况下，卫星的总能量是守恒的。由此，我们建立一个卫星轨道的能量平衡式并在嫦娥三号工程中使用这个平衡式对绕月飞行的轨道数据进行核算，最终得出近月点速度，远月点速度
以及近远月点坐标。
针对问题二，应用参数化控制求解月球探测器精确定点软着陆最优控制问题的方法。首先用约束变换技术将不等式约束进行了近似处理,而后利用若干个分段的常数去逼近最优解,再根据强化技术通过时间轴上的变换,将每一段参数的持续时间转变为一组新的参数，将最优控制问题转化为一系列参数优化问题。最后应用经典的参数优化方法求得最优控制函数的一个近似解,通过增加参数个数,重复优化得到逼近连续最优解的参数化解，并用LINGO软件最终求得，末时刻探测器质量
，燃料消耗为
，最后探测器以
的对月速度精确降落到指定登月点，并用Simulink搭建模型仿真。另外，在避障阶段，我们使用简单的区域搜索方法进行粗避障和精避障，利用MATLAB软件确定点坐标斜率大小，并画出区域搜索状况，从而避开陨石坑。
针对问题三，我们在对月球探测器向月着陆轨道的详细分析的基础上，采用协方差分析方法对其轨道初始误差及误差源造成的轨道误差进行了分析，分析了初始轨道位置误差和初始速度误差
对轨道终点的位置和速度误差的影响，及从近月点开始至进入着陆轨道为止相应的轨道位置和速度总误差的时间历程并通过问题二的答案用EXCEL画出表格。
关键词：能量守恒 软着陆 最优控制 区域搜索 误差分析
嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m。
嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段：主减速、快速调整、粗避障、精避障、减速下降、自由落体阶段。要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。
根据上述的基本要求，请你们建立数学模型解决下面的问题：
（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的 大小与方向。
（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
（3）对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
1、从飞行器在近月点下降时，月球表面是水平的。
2、不考虑空间飞行器上各点因燃料消耗而产生的位移； 3、卫星的观测图片及数据精准；
4、卫星和空间飞行器的运动是在真空中进行的； 5、卫星只受重力影响，空间飞行器除自身推力外只受重力影响； 6、在对卫星和空间飞行器进行轨道估计时，认为作用于其上的所有外力都通过
其质心；
为近地点地心距
为近月点月心距
为远地点地心距
为远月点月心距
为近地点速度
为近月点速度
为远地点速度
为远月点速度
为月球自转角速度
为月球半径
为制动火箭的比冲
为软着陆末时刻
为探测器初始质量
为制动发动机最大推力
为月球质量
为地球质量
地球卫星入轨后，处于地球的保守场中，在不计阻力损耗情况下，卫星的总能量是守恒的。由此，可以建立一个卫星轨道的能量平衡式。
卫星位置的势能为，卫星具有的动能
，则卫星的总能量为
总能量在轨道上的任何一点都是相同的。
卫星的近地点和远地点的总能量分别为
（1）
（2）
其中，GM 为地球常数，，
为近地点地心距，
为远地点地心距，
为近地点速度，
为远地点速度。
由于,有
即
（3）
由开普勒定律可知，近地点和远地点速度之比等于其地心距的反比，即
代入式（3），有
整理后，得
令有
（4）
称上式为能量平衡式，它表示在近地点动能等于其势能的倍。
若令，则有
嫦娥三号减速进入绕月飞行轨道，根据资料显示，嫦娥三号卫星进入月球引力场时，并没有形成环绕月球的轨道，嫦娥三号进入近地点高度200公里、远地点高度约38万公里的地月转移轨道。经第一次制动，直接进入了100km的圆形运行轨道。
下面按这些初始条件应用能量平衡式计算速度增量。
已知数据：
=7.3477×1022kg
=5.98×1024 kg
=1737.646km
由，则
第一次制动后达到目标轨道是(100km/100km)的圆形轨道。即
式中，为近月点月心距，
为远月点月心距。利用能量平衡式，
式中，为近月点速度，
为远月点速度。即变轨后，卫星进入了绕月的圆形轨道，速度
。
在远月点制动，目标轨道高度为（15/100） km。即
利用能量平衡式，
变轨后，近月点速度为：
远月点速度为：
由附件2知，“嫦娥”三号在软着陆过程的第一个阶段呈近抛物线运动，
抛物线的长度公示：
另外
联立方程有
利用MATLAB解方程得 。
由于能力有限，我们不考虑月球表面为弧面的影响，那么近月点就是在着陆点（19.51°W，44.12°N）上方15km为圆心，21.26km为半径的圆的点集。(事实上，当“嫦娥”三号从地球上发射时，在绕地飞行期间，速度方向就已经确定，从而在月球上的运行轨道也已经确定，但是根据给出的资料我们无法确定具体的方向以及轨道，因此将所有可能的近月点均列出来，即为下图中的圆） 如图所示：
图1近月点点集与软着陆轨道图
其中表示抛物曲线在地面的投影与着陆点所在经度的夹角。由于所给出的资料有限，我们无法确定轨道所在平面的具体位置。
若则由三角函数定理得:
则近月点的坐标为（19.60°W，44。29°N），远月点（160.4°E，44.29°S。）
若则近月点的坐标为（19.51°W，44。32°N），远月点为(160.40°E，44.32°S)。
探测器经过环月轨道的着陆方式因其具有较长的软着陆准备时间、对着陆位置的限制比较小以及减少着陆舱部分的燃料消耗等优点故而被广泛采用。该方式的关键环节就是从距离月面15 km 的近月点至月面的动力下降过程。考虑到月球自转和侧向运动,针对三维空间内精确定点软着陆问题利用参数化控制解决变推力软着陆最优控制问题。
根据资料显示软着陆分为六个阶段：主减速、快速调整、粗避障、精避障、减速下降、自由落体阶段。由于能力有限我们只讨论了避障过程，我们又将避障过程分为四个任务段: 接近段、悬停段、避障段和缓速下降段, 分别实现粗避障、高精度三维成像、精避障和着陆位置保持功能, 形成了大范围粗避障、小范围精避障和着陆位置保持的接力避障过程。
如图1 所示,定义惯性坐标系, 原点在月心,参考平面是月球赤道面,
轴指向月球赤道相对于白道的升交点,
轴指向月球自转角速度方向,
轴按右手坐标系确定。再定义月固坐标系
,
, 以月球赤道面为参考平面，
轴指向赤道面与起始子午面的交线方向，
指向月球自转角速度方向，
轴按右手坐标系确定。
为原点
图2 坐标系示意图
在探测器质心的轨道坐标系, 指向从月心到着陆器的延伸线方向,
垂直
指向运动方向,
按右手坐标系确定。制动发动机推力P的方向与探测器纵轴重合，为P与
轴正向所成夹角,Ψ为P在
平面上的投影与
轴负向所成夹角。
为
与
所成夹角,
为
在
平面上的投影与
轴正向所成夹角。
为月球自转而产生的月固坐标系相对惯性坐标系的转角,不妨假设初始时刻月固坐标系与惯性坐标系重合。
显然有轨道坐标系到惯性坐标系转换矩阵
惯性坐标系到月固定坐标系的转换矩阵为
根据牛顿第二定律，结合科氏定律整理可以得到探测器的月固坐标系中的运动方程为
其中，
，
为探测器速度矢量在月固坐标系各轴上的投影，F为发动机推力，m为探测器质量，
，
，
为该高度月球重力加速度在月固坐标系各轴上的投影，
为月球自转角速度。
因此,在月固坐标系中探测器的运动方程可表示如下：
(5)
其中
为月球引力常数,
为制动火箭的比冲,是一个常值。
取为系统状态变量,
为控制变量,则式(5) 可以简记为
按照耗燃最优的要求,取性能指标为
在实际情况下,通常没必要令探测器着陆速为零,只要能保证探测器以很小的相对速度降落到月面即可。因此,本文将软着陆的末速度要求以惩罚因子的形式加入到指标中如下式所示,主要目的是降低最优控制问题求解的复杂度,该惩罚因子可以通过反复的数值仿真运算,按经验设定。
(6)
此外,显然有约束条件
(7)
其中,
,
为预定着陆点在月固坐标系中的坐标;
为着陆点到月心距离,即月球半径。
对于含有形如这类关于状态变量在连续时间上都要满足的不等式约束最优化问题 , 至今还是最优化领域的一个难点。
显然≥ 0 等价于
（8）
但上式显然在= 0 时不可微 , 因此用如下不等式去近似上式
（9）
其中
,
是调节参数。已知当
足够小的时候 , 存在
, 使得对任何满足
的
能够令 (9) 对 (8) 达到满足要求的近似。不妨记
为用 (9) 式替换
后得到的新的约束函数。
因此本文所讨论的软着陆耗燃最优问题转化为模型2。
在系统(5) 满足约束函数。的情况下，求取适当的控制变量
使指标函数(6) 达到最小。
假定初始时刻为0 ,终端时刻为待定参数。选取满足
的序列
和三组参数
， 构造形如
的参数化分段常数控制器。其中
用控制器替换系统(4) 中的
, 则模型2转变为模型3。
寻找三组参数来最小化指标函数(5) ,并且满足约束函数
。
显然,对于每个给定的,这都是一个有限维的参数优化问题。已知当
时,模型3的最优解收敛于模型2的最优解。不过前辈已经证明了在数值计算中,求解模型3的参数梯度时难度很大甚至求不出真实解,因而本文引入强化技术来解决这一问题。
从s ∈[0 ,1 ] 到t ∈[0 ,构造如下
上式中，
，序列
为[0,1]区间上预先给定的分段点，并且满足
。
将上式两边对s求导可得
其中
不妨令
则得到如下增广系统
即
(10)
其中,
,
与
,
,
分别为
,
,
与
,
,
经过变换后的形式，
，
，
，
，
，
,
。
指标函数变为
(11)
约束条件变为
(12)
其中。
由于仅仅已知探测器在软着陆起始点到月心的距离和探测器的初始速度以及初始质量
， 而软着陆起始点另两个空间位置信息
角与β角的初始值
与
未知，因而令
与
为系统待定参数，则系统初始状况可以表示为
(13)
那么模型3转化为如下模型4。
在系统(10)满足约束并且初始条件如式(13)的情况下,求取适当的控制变量使指标函数(11)达到最小。
再由
可知,模型4将最初的探月飞行器软着陆最优控制问题转化成了优化静态控制参数，
，
和
以及系统参数
，
，
的问题，利用经典的参数优化算法即可求出登月飞行器的软着陆最优控制的一组逼近解和软着陆最优初值点位置以及终端时刻。利用此算法，增加时间的分段点个数可以重新优化，经过多次优化后即可得到满意精度的参数化解。
此外,假如令系统(5)中的推力F 为已知的恒定推力，令控制变量,则本文问题变为恒定推力下软着陆最优控制问题，依然可以利用本文方法解决,而依据极大值原理结合传统的打靶法则只能解决恒定推力的情况。
设探测器初始质量，制动发动机最大推力
，比冲
；初始速度
，
，
；月球自转角速度
；月球引力常数
；近月点距月心距离
，月球半径
。着陆点位置为西经19.51°，北纬44.12°,。
利用LINGO软件，通过计算机运算，令，
，
即可得到符合精度的最优解，最终利用本文的参数化控制得到软着陆末时刻
，末时刻探测器质量0.846t，燃料消耗为1.554t ，最后探测器以3. 615 m/ s 的对月速度精确降落到指定登月点。
根据粗避障任务的要求，飞行轨迹要保证：1)成像敏感器能够持续观测预定着陆区范围； 2) 着陆器无机动就能到达预定着陆区。考虑到7500 N 主发动机羽流带来的不可见区域为半锥角约25°的锥体，而成像敏感器视场为30°，为了避免主发动机羽流对成像敏感器的影响，且使成像敏感器的视线距离尽可能短，取成像敏感器视线偏置40°角。为了保证在接近段成像敏感器视场能够观测到着陆区，确定采用下降轨迹接近与水平面夹角45°的直线下降方式逐步接近着陆区.。
悬停段主要任务就是对月面进行高精度三维成像， 精确检测着陆区域的障碍，确定安全着陆点。因此， 着陆器需要保持悬停状态，速度和姿态角速度接近0值，姿态和位置不变。根据激光三维成像敏感器的工作范围限制和观测足够大着陆区的要求，选取高度在离着陆点上方100m处悬停，悬停状态下利用激光三维成像敏感器精确观测着陆区，并处理三维图像数据，确定安全着陆点，转入精避障段。
图3 距月面2400m处的数字高程图
在这个粗避障过程中我们用到简单的区域搜索模型：在粗避障阶段在检测到的某一点X为中心，取出的矩阵
.
求出每一点的斜率。控制
。
若则记为平坦，否则记为不平坦。以此可得到一个矩阵Y。由此我们利用MATLAB编程得到如下图所示：
图4 粗避障区域搜索图
为了保证最后的避障落点精度和节省推进剂，着陆器精确避障和下降同时进行。并悬停在距离月面100m处，对着陆点附近区域100m范围内拍摄图像，并获得三维数字高程图。 根据确定的安全着陆点，从约100 m高处斜向下降到着陆点上方30 m，相对月面下降速度到1.5 m/s，水平速度接近0。
图5 距离月面100m 处的数字高程图
与粗避障相似，同样用区域搜索模型，并利用MATLAB编程得出如下搜索图：
图6 精避障区域搜索图
为了保证着陆月面的速度和姿态控制精度, 缓速下降段要以较小的设定速度匀速垂直下降, 消除水平速度和加速度, 保持着陆器水平位置, 直到收到关机敏感器测量信号或加速度测量大于预设置值，就关闭主发动机.。即实现在距离月面4m处相对月面静止，使嫦娥三号自由落体到精确有落月点。考虑到推进剂的消耗和导航位置漂移，选择下降速度为-2 m/s(选定月心到飞行器的方向为正方向)。
引起轨道误差的误差源主要是速度方向误差和导航误差，包括位置误差和速度误差。其中：位置误差：，分别为在月心惯性坐标系中
轴、
轴、
轴的分量。速度误差：
分别是在月心惯性坐标系
轴、
轴、
轴的分量。向月飞行着陆轨道的初始轨道位置和速度误差由嫦娥三号的发射在轨道上的精度决定，若探测器在飞行途中进行轨道修正，则经过轨道时的轨道位置误差将由导航误差决定，速度误差将由姿态误差和制导误差决定。上述误差决定了轨道误差协方差分析的计算初始条件，表1 给出了在不进行中途轨道修正情况下，在月心惯性坐标系里，初始轨道位置（即近月点位置）误差和初始速度误差对轨道终点的位置和速度误差的影响。
表1 初始轨道位置和速度误差对轨道终点误差的影响
6.1 模型评价
在问题一中由于地球卫星入轨后，处于地球的保守场中，在不计阻力损耗情况下，卫星的总能量是守恒的。由此，我们建立的一个卫星轨道的能量平衡式。在“嫦娥”三号工程中，使用这个平衡式对绕月飞行的轨道数据进行核算是比较准确的，而且物理概念也很清晰。
问题二中的模型相比较传统的打靶法而言,我们的方法无需猜测毫无物理意义的共轭变量初值,既可以解决变推力软着陆问题,又适用于恒定推力的情况,同时优化速度更快。
但在结果计算中，部分数据理想化(如经纬度计算中忽略了经线位置不同时度数的计算也不相同以及令发动力始终最大为)，致使结果不够精准。
第一，由于在上述模型中没有考虑到空间飞行器上各点因燃料消耗而产生位移的情况。可以建立燃料消耗与位移的关系模型来确定飞行器的位移变化。 第二，从飞行器在近月点下降时，可以准确的考虑月球表面的弧度，这样可以更精确的确定近月点与远月点的坐标。 第三，在上述的基础上，可以用基于蚁群算法和基于广义乘子法等对月球软着陆轨道的优化设计进行比较，得出更精确的着陆轨道。
[1]郗晓宁，王威。近地航天器轨道基础[M]。长沙: 国防科技大学出版社，2003。[2]叶培建，孙泽洲。嫦娥工程技术手册（卫星分册）[R]。北京：国防科工委月
球探测工程中心，2007。[3]朱建丰，徐世杰。基于自适应模拟退火遗传算法的月球软着陆轨道优化[J] 航空学报，2007，28。[4]段佳佳，徐世杰，朱建丰。基于蚁群算法的月球软着陆轨迹优[J]。宇航学报，
2008，29。[5]屠善澄。卫星姿态动力学与控制。北京：宇航出版社，1998。168–175。[6]王大轶，李铁寿，马兴瑞。月球最优软着陆两点边值问题的数值解[J]。航天
控制，2000，3：44-49。[7]赵福安，王长钰。非线性规划问题的灵敏度分析和稳定性分析 http://www.docin.com/p-209433796.html 2014年09月14
附录
Matlab运行代码：
Untitled03.m
z=imread('F:\qqfile\2014\A\??3 ?2400m???????.tif')
x=1:2300;
[X,Y]=meshgrid(x);
Z=double(z);
mesh(X,Y,Z)
Untitled4.m
z=imread('F:\qqfile\2014\A\??4 ???100m???????.tif')
x=1:1000;
[X,Y]=meshgrid(x);
Z=double(z);
mesh(X,Y,Z)
ceshi.m
X=imread('F:\qqfile\2014\A\??3 ?2400m???????.tif');
R=10; %????
k=tan(pi/6); %????,????????pi/6
[ Xe ] = roughAvoid( X, R, k ); %?????,????Xe
imwrite(Xe,'Xe.bmp','bmp'); %?Xe?????Xe.bmp??
Ceshi1.m
X=imread('F:\qqfile\2014\A\??4 ???100m???????.tif');
R=4; %????
k=tan(pi/6); %????,????????pi/6
[ Xe1 ] = roughAvoid( X, R, k ); %?????,????Xe
imwrite(Xe1,'Xe1.bmp','bmp'); %?Xe?????Xe.bmp??
isFlat.m
function [ C ] = isFlat( Y, k )
% ??Y??????,??????1,?????,??0
%
[m,n]=size(Y);
C=ones(size(Y));
Yl=Y(:,1);
Yr=flipud(Y(:,m));
for i=1:m
j=m-i+1;
kc=(Y(j,m)-Y(i,1))/sqrt((j-i)^2+m^2);
if abs(kc)>k
C=zeros(size(Y));
return;
end
end
for i=1:n
j=n-i+1;
kc=(Y(n,j)-Y(1,i))/sqrt((j-i)^2+n^2);
if abs(kc)>k
C=zeros(size(Y));
return;
end
end
end
roughAvoid.m
function [ Xe ] = roughAvoid( X, R, k )
% ???,???X???????2R+1?????????
% ??X?????(m.,n)?????(2R+1)*(2R+1)??X(m-R:m+R,n-R:n+R)
% ???
% X??
% R????????
% k??????
[mX, nX]=size(X);
Xe=ones(mX, nX);
for m=R+1:R:mX-R %????R??????
for n=R+1:R:nX-R %????R??????
[ Xe(m-R:m+R,n-R:n+R) ] = isFlat( X(m-R:m+R,n-R:n+R), k );
end
end
end
roughAvoid1.m
function [ Xe1 ] = roughAvoid( X, R, k )
% ???,???X???????2R+1?????????
% ??X?????(m.,n)?????(2R+1)*(2R+1)??X(m-R:m+R,n-R:n+R)
% ???
% X??
% R????????
% k??????
[mX, nX]=size(X);
Xe1=ones(mX, nX);
for m=R+1:R:mX-R %????R??????
for n=R+1:R:nX-R %????R??????
[ Xe1(m-R:m+R,n-R:n+R) ] = isFlat( X(m-R:m+R,n-R:n+R), k );
end
end
end
Lingo运行代码：
sets: variables/1..20/:afa,t,v,r,m,theita,gama,afac,tc,vc,rc,mc,theitac,gamac,fHv,fHr,fHm,fHtheita,fHgama;endsetsdata:R0=1752646; g0=1.63;ge=9.8; Isp=300;m00=2400; VCC=1692.00;n=0.01365; afa0=3.615;t0=0;v0=1.049;r00=1.0094;m0=1;theita0=0;gama0=-0.0436;enddata@for(variables(i):@free(afa(i));@free(t(i));@free(v(i));@free(r(i));@free(m(i));@free(theita(i));@free(gama(i));@free(afac(i));@free(tc(i)); @free(vc(i));@free(rc(i));@free(mc(i));@free(theitac(i));@free(gamac(i)); @free(fHv(i));@free(fHr(i));@free(fHm(i));@free(fHtheita(i));@free(fHgama(i)); @free(gama0););[obj]min=-m(20);!ç؛؟و€§ن¸چç‰ه¼ڈç؛¦وں;@for(variables(i):@bnd(1.57,afa(i),4.71); @bnd(0,t(i),3000); @bnd(1.57,afac(i),4.71); @bnd(0,tc(i),3000); r(i)>=1;);!ç؛؟و€§ç‰ه¼ڈç؛¦وں;r(20)=1;v(20)=10.655/VCC; !éç؛؟و€§ç؛¦وں;!i=1و—¶;@for(variables(i)|i#eq#1:vc(i)=(v0+v(i))/2+n/8*(-@sin(gama0)/r00^2+t0/m00/g0*@cos(afa0)/m0+@sin(gama(i))/r(i)^2-t(i)/m00/g0*@cos(afa(i))/m(i)); rc(i)=(r00+r(i))/2+n/8*(v0*@sin(gama0)-v(i)*@sin(gama(i))); mc(i)=(m0+m(i))/2+n/8*(-@sqrt(R0/g0)*t0*r00^2/Isp/m00/g0+@sqrt(R0/g0)*t(i)*r(i)^2/Isp/m00/g0); theitac(i)=(theita0+theita(i))/2+n/8*(v0*@cos(gama0)/r00-v(i)*@cos(gama(i))/r(i)); gamac(i)=(gama0+gama(i))/2+n/8*((v0^2/r00-1/r00^2)*@cos(gama0)/v0+t0/m00/g0*@sin(afa0)/m0/v0-(v(i)^2/r(i)-1/r(i)^2)*@cos(gama(i))/v(i)-t(i)/m00/g0*@sin(afa(i))/m(i)/v(i)); fHv(i)=3/2/n*(v(i)-v0)-(-@sin(gama0)/r00^2+t0/m00/g0*@cos(afa0)/m0-@sin(gama(i))/r(i)^2+t(i)/m00/g0*@cos(afa(i))/m(i))/4; fHr(i)=3/2/n*(r(i)-r00)-(v0*@sin(gama0)+v(i)*@sin(gama(i)))/4; fHm(i)=3/2/n*(m(i)-m0)-(-@sqrt(R0/g0)*t0*r00^2/Isp/m00/g0-@sqrt(R0/g0)*t(i)*r(i)^2/Isp/m00/g0)/4; fHtheita(i)=3/2/n*(theita(i)-theita0)-(v0*@cos(gama0)/r00+v(i)*@cos(gama(i))/r(i))/4; fHgama(i)=3/2/n*(gama(i)-gama0)-((v0^2/r00-1/r00^2)*@cos(gama0)/v0+t0/m00/g0*@sin(afa0)/m0/v0+(v(i)^2/r(i)-1/r(i)^2)*@cos(gama(i))/v(i)+t(i)/m00/g0*@sin(afa(i))/m(i)/v(i))/4;);! i=2هˆ°20و—¶;!هگ„点ن¸é—´هڈکé‡ڈه€¼;@for(variables(i)|i#ge#2:vc(i)=(v(i-1)+v(i))/2+n/8*(-@sin(gama(i-1))/r(i-1)^2+t(i-1)/m00/g0*@cos(afa(i-1))/m(i-1)+@sin(gama(i))/r(i)^2-t(i)/m00/g0*@cos(afa(i))/m(i)); rc(i)=(r(i-1)+r(i))/2+n/8*(v(i-1)*@sin(gama0)-v(i)*@sin(gama(i))); mc(i)=(m(i-1)+m(i))/2+n/8*(-@sqrt(R0/g0)*t(i-1)*r(i-1)^2/Isp/m00/g0+@sqrt(R0/g0)*t(i)*r(i)^2/Isp/m00/g0); theitac(i)=(theita(i-1)+theita(i))/2+n/8*(v(i-1)*@cos(gama(i-1))/r(i-1)-v(i)*@cos(gama(i))/r(i)); gamac(i)=(gama(i-1)+gama(i))/2+n/8*((v(i-1)^2/r(i-1)-1/r(i-1)^2)*@cos(gama(i-1))/v(i-1)+t(i-1)/m00/g0*@sin(afa(i-1))/m(i-1)/v(i-1)-(v(i)^2/r(i)-1/r(i)^2)*@cos(gama(i))/v(i)-t(i)/m00/g0*@sin(afa(i))/m(i)/v(i)); fHv(i)=3/2/n*(v(i)-v(i-1))-(-@sin(gama(i-1))/r(i-1)^2+t(i-1)/m00/g0*@cos(afa(i-1))/m(i-1)-@sin(gama(i))/r(i)^2+t(i)/m00/g0*@cos(afa(i))/m(i))/4; fHr(i)=3/2/n*(r(i)-r(i-1))-(v(i-1)*@sin(gama0)+v(i)*@sin(gama(i)))/4; fHm(i)=3/2/n*(m(i)-m(i-1))-(-@sqrt(R0/g0)*t(i-1)*r(i-1)^2/Isp/m00/g0-@sqrt(R0/g0)*t(i)*r(i)^2/Isp/m00/g0)/4; fHtheita(i)=3/2/n*(theita(i)-theita(i-1))-(v(i-1)*@cos(gama(i-1))/r(i-1)+v(i)*@cos(gama(i))/r(i))/4; fHgama(i)=3/2/n*(gama(i)-gama(i-1))-((v(i-1)^2/r(i-1)-1/r(i-1)^2)*@cos(gama(i-1))/v(i-1)+t(i-1)/m00/g0*@sin(afa(i-1))/m(i-1)/v(i-1)+(v(i)^2/r(i)-1/r(i)^2)*@cos(gama(i))/v(i)+t(i)/m00/g0*@sin(afa(i))/m(i)/v(i))/4;);@for(variables(i): -@sin(gamac(i))/rc(i)^2+tc(i)/m00/g0*@cos(afac(i))/mc(i)=fHv(i); vc(i)*@sin(gamac(i))=fHr(i); -@sqrt(R0/g0)*tc(i)*rc(i)^2/Isp/m00/g0=fHm(i); vc(i)*@cos(gamac(i))/rc(i)=fHtheita(i); (vc(i)^2/rc(i)-1/rc(i)^2)*@cos(gamac(i))/vc(i)+tc(i)/m00/g0*@sin(afac(i))/mc(i)/vc(i)=fHgama(i); );
本文来源：https://www.2haoxitong.net/k/doc/70fe34def01dc281e53af0b2.html
文档为doc格式