嫦娥三号软着陆最优控制建模与Python求解

简介一份2014年全国大学数学建模竞赛A题完整论文主题为嫦娥三号软着陆轨道设计与控制策略适合参加数模竞赛的学生、航天工程爱好者以及需要完成相关毕业设计的读者参考。论文围绕着陆准备轨道近月点与远月点位置、速度方向求解以及六阶段软着陆轨道和控制策略优化展开。采用逆向推理与微元分析方法计算出探测器水平位移514.8km确定近月点19.51W,27.08N上空15km、远月点19.51E,152.92S上空100km的位置通过模拟退火算法得到最小燃料消耗468.25kg并绘制安全区域与软着陆轨道图像最后进行误差分析和敏感性分析。资源为1个Word文档.doc压缩包体积771KB论文结构完整包含摘要、问题重述、模型假设、模型建立与求解等部分可直接作为数学建模论文写作的参考模板。目前已有87人学习对掌握嫦娥三号软着陆轨道设计方法和优化技巧很有帮助。1. 嫦娥三号软着陆问题的数学本质2014年全国大学生数学建模竞赛A题把“嫦娥三号软着陆”拆成两个动作轨道设计与控制策略。前者离线算出一条从15公里停泊轨道到月面落点的可行轨迹后者解决实际飞行中发动机推力大小与方向如何随时间变化。抽象成数学问题这是一个带终端速度约束、推力幅值硬约束与燃料积分约束的非线性最优控制问题。很多参赛队拿到题目后先去搜嫦娥三号任务本身的参数却忽略了题目给出的初始高度、落点坐标和发动机参数才是建模的输入。这个问题的核心难点不在天体力学而在如何处理“不可关机的发动机”和“终端速度归零”这一对矛盾。本文按建模、求解、制导、写作四步展开提供可直接改参数运行的Python代码最后补充评阅视角下最容易被扣分的两个细节。2. 着陆动力学模型与坐标系的选取2.1 二维垂直平面模型的几何设定嫦娥三号软着陆过程从距月面约15公里的停泊点开始初始速度方向接近水平量级约为1692米每秒。若完整使用三维方程状态量包含三个位置、三个速度和质量共七个变量对于竞赛论文来说公式量和调试成本都会翻倍。常见做法是将问题限制在“包含月心、初始速度方向和目标落点”的二维垂直平面内只保留水平位移x和高度z两个位置量。这个二维近似的成立条件是整个下降段的航程不超过数十公里落点经度与纬度变化对重力加速度方向的影响小于1%。在这个尺度下月球表面可以视为平面重力加速度取常数1.62米每二次方秒。这样做不会带来实质性的模型失真却能让控制量从三维推力方向角简化为一个俯仰角θ。2.2 状态方程与质量消耗选择以落点为原点的地平直角坐标系x轴为水平方向、z轴垂直向上。状态向量定义为X [x, z, vx, vz, m]其中m为着陆器当前质量。推力幅值F与俯仰角θ作为两个控制量θ表示推力方向与水平面的夹角单位弧度。运动方程为水平方向加速度由推力水平分量产生垂直方向加速度由推力垂直分量与月球重力共同决定质量变化率由比冲和推力幅值决定。用Python封装这组常微分方程方便后续积分器复用import numpy as np G_MOON 1.62 # 月球表面重力加速度m/s² ISP 3000.0 # 发动机比冲秒 VE ISP * 9.8 # 等效排气速度m/s def landing_ode(t, X, F, theta): 动力下降段状态方程。 X [x, z, vx, vz, mass] theta 为推力与水平面夹角向上为正单位弧度。 返回 dX/dt。 x, z, vx, vz, mass X c, s np.cos(theta), np.sin(theta) ax F * c / mass az F * s / mass - G_MOON dm -F / VE return np.array([vx, vz, ax, az, dm])参数说明推力方向角θ接近90度表示推力基本垂直向上用于悬停或抵消重力θ较小时推力沿水平方向作用用于削减水平速度。比冲Isp越高同样推力下燃料消耗越慢。若题目给定的是推力与比冲的范围这里可以取中值作为标称参数后续做敏感性分析时再扰动这两个值。2.3 控制约束与边界条件表建模型时最容易忽略的是发动机“不可关机”这一工程约束。液体发动机或单组元推进器存在燃烧稳定下限推力不能低于某个百分比一般取最大推力的20%作为F_min。这个下限会让最优解从“开关控制”变成“常值推力弧段”对轨迹形态影响很大。约束类型数学表达物理含义推力幅值下限F ≥ F_min 0.2 × F_max发动机不可关机推力幅值上限F ≤ F_max结构强度限制推力方向归一化cos²θ sin²θ 1俯仰角定义初始状态x0, z15000, vx1692, vz0停泊点终端位置x(t_f)0, z(t_f)0到达目标落点终端速度vx(t_f)≈0, vz(t_f)≈0软着陆判定路径约束z(t) ≥ 0全程成立不撞月面F_min不为零这一点直接决定了求解方法的选择。若允许推力归零最优控制问题的开关函数会出现奇异弧段而加上F_min后最常见的工程解是主制动段全程满推力末段切换到小推力精确调整。后文的三段式策略正是围绕这个约束展开的。2.4 三段式着陆过程的划分工程上把整个软着陆过程分成主制动段、接近段和垂直下降段。主制动段从15公里高度开始占全部燃料消耗的八成以上发动机保持近最大推力推力方向与速度方向相反并略偏上目的是把接近1.7公里每秒的水平速度降到几十米每秒量级。接近段从高度约3公里开始此时水平速度已小发动机推力下调至接近悬停值推力方向逐渐转向垂直把下降速率稳定在每秒2米左右。垂直下降段从高度约300米开始推力基本指向上方先悬停后缓慢下降直至触地关机。这三段的切换条件在第四章节用闭环状态机实现。每段内部的制导规律不同但共用同一个landing_ode动力学函数积分器不需要区分阶段。3. 燃料最优控制模型与直接打靶法求解3.1 从性能指标到最优控制问题主制动段的燃料消耗正比于推力对时间的积分。由于这一段的推力恒定在最大值附近最小化燃料等价于最小化主制动段时间。但整个软着陆过程包含末段的小推力调整不能只用最短时间做全局指标因此写成更一般的形式目标是极小化J等于从零到终端时刻T对推力F(t)的积分。约束包括状态方程、边界条件和推力幅值范围。这个最优控制问题的哈密顿函数中推力方向只通过速度伴随变量进入因此最优推力方向应该与速度伴随向量反平行。这条结论说明一个问题推力方向不能凭直觉设定为“与速度方向相反”因为伴随变量和状态速度在终端并不平行。对应到工程表达就是主制动段横向减速和垂直减速之间存在一个最优的比例分配这个比例随时间变化而不是固定的反推角。使用间接法求这个最优比例需要猜测伴随变量初值对新手不友好竞赛场景下更稳妥的是直接法。3.2 直接打靶法及其数值特性直接打靶法的做法是把整个时间区间等分成N段每段内的控制量设为常数然后将终端误差作为优化目标用数值优化算法迭代调整控制序列。相比间接法不存在伴随变量初值猜测的问题比直接配点法少一整套非线性规划求解器的配置是竞赛中最容易复现的方案。优化变量包括每段的俯仰角θ₁到θ_N以及总时长T。主制动段推力固定为F_max。目标函数由终端位置误差、终端速度误差和相邻控制量的光滑惩罚组成。振动所在若不使用光滑惩罚控制序列会出现逐段抖振若加大光滑权重终端精度又会受损。求解方法精度实现难度对初值的敏感度竞赛推荐度间接共轭梯度极高高非常敏感不推荐直接打靶法中高低中等初值需物理合理推荐直接配点法高高稳健有余力时选用直接打靶法的另一个优势是它天然与仿真代码共用同一套动力学函数写完优化器之后把最优控制序列交给闭环仿真去跑不需要重新推导模型。3.3 直接打靶法参考代码以下代码先将时间离散成40段每段常值俯仰角再调用scipy的BFGS优化器寻找终端误差最小的控制序列。from scipy.optimize import minimize from numpy import deg2rad, diff, hypot N_SEG 40 F_MAX 7500.0 X0 np.array([0.0, 15000.0, 1692.0, 0.0, 1200.0]) def rk4_step(f, t, X, h, *args): k1 f(t, X, *args) k2 f(t h/2, X h/2 * k1, *args) k3 f(t h/2, X h/2 * k2, *args) k4 f(t h, X h * k3, *args) return X h * (k1 2*k2 2*k3 k4) / 6.0 def simulate(theta_seq, T_total): dt T_total / N_SEG X X0.copy() t 0.0 for k in range(N_SEG): X rk4_step(landing_ode, t, X, dt, F_MAX, theta_seq[k]) t dt return X def objective(p): theta_seq p[:N_SEG] T_total p[N_SEG] xf, zf, vxf, vzf, mf simulate(theta_seq, T_total) # 特征尺度归一化位置、速度、光滑项都在0~1量级 pos_err (xf / 3000.0) ** 2 (zf / 100.0) ** 2 vel_err (vxf / 50.0) ** 2 (vzf / 20.0) ** 2 smooth 0.05 * np.sum(diff(theta_seq) ** 2) return pos_err vel_err smooth theta0 np.linspace(deg2rad(65), deg2rad(88), N_SEG) p0 np.r_[theta0, 600.0] res minimize(objective, p0, methodBFGS, options{maxiter: 300, ftol: 1e-8}) print(最优总时长:, res.x[N_SEG], 秒)代码说明目标函数中的pos_err和vel_err分别把位置误差除以3000米和100米、速度误差除以50米每秒和20米每秒目的是让误差量纲统一到同一数量级。如果不做这一步位置误差以米计、速度误差以米每秒计数值优化会将全部权重偏向位置项得到一个落点准确但终端速度极大的假解。光滑项系数0.05在工程调试中需要根据仿真结果调整系数过大会使轨迹偏柔但落点偏差变大。3.4 优化失败时的诊断方法直接打靶法最常见的两个失败场景一是总时长T被优化到极大值轨迹出现明显“飘”的现象原因是目标函数没有加入路径约束z(t)≥0的惩罚着陆器在优化的中间迭代中穿过月面以下终端又飞回零点。解决办法是给目标函数增加一项max(0, -z_min)^2的罚函数z_min在仿真中记录即可。二是推力方向初始猜测过于小例如从30度开始优化水平推力占主导着陆器在水平方向漂移几公里位置误差超过特征归一化尺度导致梯度失效。经验做法是θ初值范围内65度和88度之间线性插值初始阶段偏水平减速末段偏垂直下降与物理直觉一致。若终端速度还降不下来把目标函数中速度项的特征尺度调小逼迫优化器优先保证速度为零。4. 闭环制导控制策略与仿真参数调优4.1 重力转弯与位置修正结合的制导律开环最优轨迹求出来之后还需要一套闭环控制律应对模型误差和干扰。主制动段常见做法是“重力转弯制导”即推力方向跟踪速度反方向同时叠加一个指向目标点的位置修正项。修正权重不需要很大8%左右就能避免落点漂移太大则会引起轨迹振荡。垂直方向的控制则用比例反馈。接近段的目标是把垂直速度稳在每秒2米的下落水平速度压到接近零垂直下降段把目标改为每秒1米到月面附近再缓缓逼近零。反馈增益的选择以临界阻尼为目标即用0.2到0.5之间的速度增益配合重力加速度前馈避免着陆器在接近月面时出现明显的上下抖动。以下是主制动段与接近段的闭环控制函数供仿真循环直接调用def closed_loop_step(state, phase, dt): x, z, vx, vz, m state if phase brake: vnorm np.hypot(vx, vz) 1e-6 u_speed -np.array([vx, vz]) / vnorm u_pos np.array([-x / 5000.0, 0.0]) u u_speed 0.08 * u_pos u u / np.linalg.norm(u) theta np.arctan2(u[1], u[0]) F F_MAX elif phase approach: # 水平速度反馈归零垂直速度收敛到 -2 m/s Fx m * 0.02 * (0 - vx) Fz m * (G_MOON 0.25 * (-2.0 - vz)) F np.clip(np.hypot(Fx, Fz), 1500.0, 4500.0) theta np.arctan2(Fz, Fx) else: # 垂直下降段悬停 缓慢下降 Fz m * (G_MOON 0.35 * (-1.0 - vz)) F np.clip(Fz, 1200.0, 2500.0) theta np.pi / 2.0 X_next rk4_step(landing_ode, 0.0, state, dt, F, theta) return X_next逻辑说明制动段的u_pos取落点横向偏差除以5000米进行归一化修正权重0.08经过调参得到若落点偏差超过1公里这个修正项会把推力方向额外偏转约9度可有效拉回轨道。接近段的Fx、Fz不再直接计算俯仰角而是将需要的加速度换算成力再用反正切求出角度这种写法让物理意义更清楚也方便限制推力上限。4.2 阶段切换条件与迟滞阶段切换不能只用单一高度阈值否则仿真在阈值附近会出现频繁抖动。主制动段切换至接近段的条件是高度低于3000米且水平速度小于80米每秒接近段切换至垂直下降段的条件是高度低于300米且水平速度小于5米每秒。加入速度条件后切换变得平稳。实际仿真中还需要加入迟滞区间。例如进入接近段后若高度回到3200米以上需要退回主制动段进入垂直下降段后若垂直速度出现正值说明着陆器在向上飘不能立即切回而是继续垂直控制把上升速度压回负值。这个迟滞逻辑在状态机里用两个独立布尔变量记录即可。4.3 关键参数表与敏感性试验设计仿真参数直接决定轨迹形态这里列出主制动段、接近段和垂直下降段的标称取值。切换高度与增益并非固定的需要根据打靶法得到的开环轨迹做微调。参数主制动段接近段垂直下降段高度范围15km~3km3km~300m300m~0m推力范围5500~7500N1500~4500N1200~2500N俯仰角范围55°~88°85°~89°接近90°水平速度反馈增益0.020.020垂直速度反馈增益无0.250.35仿真完成后必须做参数敏感性分析。常见做法是把重力加速度从1.62改到1.55和1.69把比冲从3000改到2700和3300分别重新跑一次同一闭环策略。观察落点偏差和终端速度是否仍满足软着陆指标。若重力增加时终端垂直速度超过2米每秒则需要把垂直反馈增益调高若比冲下降导致燃料不足则需要向后移动接近段的切换高度。5. 把模型写成能评阅的论文三张图与两个坑参赛队伍提交的Word论文是评阅的唯一依据评阅人通常先看图、再扫公式、最后查结论。所以至少准备三张图一张是高度与水平距离的轨迹图标注停泊点、主制动段终点和落点一张是推力幅值与俯仰角随时间变化的控制指令图这张图必须能看出三段式分段和控制策略的差异第三张是敏感性分析图展示初始速度误差正负5%或推力误差正负2%时的终端落点散布。三张图分别对应轨道设计、控制策略、鲁棒性评价三个得分点。第一个坑是没有把“软着陆”工程约束转成数学约束。评阅中最常见的问题是终端速度写为零但没有提F_min下限或者路径约束z(t)≥0完全缺失。正确的做法是在模型部分明确写出带不等式约束的最优控制问题并说明数值求解时把不等式以惩罚项加入目标函数。论文里给出惩罚项的具体表达式评阅人会基于写作逻辑判断是否掌握了处理方法而不是只看仿真截图是否好看。第二个坑是不做尺度归一化竞赛论文里同时出现千米、米每秒和千克时。位置误差可能达到千米量级、速度误差只有米每秒量级两者简单加和会让优化偏向位置精度。论文中应写明“位置误差除以10公里、速度误差除以1公里每秒、质量除以初始质量”用一句话说明所有误差量统一定义在0到1区间。这个细节代表你是否真正调通过数值优化器评阅人阅读模型部分时会重点留意。最后一个进阶技巧是在仿真结果中放两组对比曲线标称模型与模型参数摄动正负10后的结果。题目给定的物理常数必然存在不确定性评阅人期望看到控制策略在偏差存在时仍能保守收敛。这个检验不需要额外建模同一套闭环控制函数更换G_MOON和ISP运行即可然后把两组结果的终端状态差做成数据表放在正文。一个数据表比单独画一条完美的标称着陆曲线更有说服力因为它证明你的策略不是对某组特定参数调出来的。本文还有配套的精品资源点击获取