http://www.ma-xy.com
目录
第一章 嫦娥 3 号 1
1.1 题目要求 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 附件 1：问题的背景与参考资料； . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2 附件 2：嫦娥三号着陆过程的六个阶段及其状态要求； . . . . . . . . . . . 4
1.1.3 附件 3：距月面 2400m 处的数字高程图； . . . . . . . . . . . . . . . . . . 6
1.1.4 附件 4：距月面 100m 处的数字高程图。 . . . . . . . . . . . . . . . . . . . 6
1.2
嫦娥
3
号软着陆的优化控制
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.1 基础知识准备 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.2 模型假设 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.2.3 符号说明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.2.4 问题一的分析与求解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.2.5 问题二的分析与求解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1.2.6 问题三的分析与求解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
1
http://www.ma-xy.com
第一章 嫦娥 3 号
1.1 题目要求
A 题：嫦娥三号软着陆轨道设计与控制策略
¬
嫦娥三号于 2013 年 12 月 2 日 1 时 30 分成功发射，12 月 6 日抵达月球轨道。嫦娥三号在
着陆准备轨道上的运行质量为 2.4t，其安装在下部的主减速发动机能够产生 1500N 到 7500N 的
可调节推力，其比冲 (即单位质量的推进剂产生的推力) 为 2940m/s，可以满足调整速度的控制
要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发
动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为 19.51W ，44.12N，海拔为
−2641m(见附件 1)。
嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着
陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点 15km，远月点
100km 的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为 6 个阶段 (见附件
2)，要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。
根据上述的基本要求，请你们建立数学模型解决下面的问题：
(1) 确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。
(2) 确定嫦娥三号的着陆轨道和在 6 个阶段的最优控制策略。
(3) 对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
1.1.1 附件 1：问题的背景与参考资料；
1. 中新网 12 月 12 日电 (记者姚培硕)
根据计划
，嫦娥三号将在北京时间 12 月 14 号在月球表面实施软着陆。嫦娥三号如何实现
软着陆以及能否成功成为外界关注的焦点。目前，全球仅有美国、前苏联成功实施了 13 次无人
月球表面软着陆。
北京时间 12 月 10 日晚，嫦娥三号已经成功降轨进入预定的月面着陆准备轨道，这是嫦娥
三号“落月”前最后一次轨道调整。在实施软着陆之前，嫦娥三号还将在这条近月点高度约 15
公里、远月点高度约 100 公里的椭圆轨道上继续飞行。期间，将稳定飞行姿态，对着陆敏感器、
¬
2014 高教社杯全国大学生数学建模竞赛题目
http://www.chinanews.com/mil/2013/12-12/5608941.shtml
1
http://www.ma-xy.com
1.1 题目要求 第一章 嫦娥 3 号
着陆数据等再次确认，并对软着陆的起始高度、速度、时间点做最后准备。
“发射、近月制动、变轨和月面降落比较起来，后者更为关键。这对我们来说是一个全新的，
也是一个最重要的考验。”中国探月工程总设计师吴伟仁表示。
嫦娥三号着陆地点选在较为平坦的虹湾区。但由于月球地形的不确定性，最终“落月”地点
的选择仍存在一定难度。据悉，嫦娥三号将在近月点 15 公里处以抛物线下降，相对速度从每秒
1.7 公里逐渐降为零。整个过程大概需要十几分钟的时间。探测器系统副总指挥谭梅将其称为“黑
色 750 秒”。
由于月球上没有大气，嫦娥三号无法依靠降落伞着陆，只能靠变推力发动机，才能完成中途
修正、近月制动、动力下降、悬停段等软着陆任务。据了解，嫦娥三号主发动机是目前中国航天
器上最大推力的发动机，能够产生从 1500 牛到 7500 牛的可调节推力，进而对嫦娥三号实现精
准控制。
在整个“落月”过程中，“动力下降”被业内形容为最惊心动魄的环节。在这个阶段，嫦娥三
号要完全依靠自主导航控制，完成降低高度、确定着陆点、实施软着陆等一系列关键动作，人工
干预的可能性几乎为零。“在这个时间段内测控都跟不上了，判断然后上去执行根本来不及，只
能事先把程序都设定好。”谭梅表示。
在距月面 100 米处时，嫦娥三号要进行短暂的悬停，扫描月面地形，避开障碍物，寻找着陆
点。“如果下面有个大坑，需要挪个地方，它就会自己平移，等照相机告诉它地面平了，才会降
落”。中国绕月探测工程首任首席科学家、中国科学院院士欧阳自远介绍。
之后，嫦娥三号在反推火箭的作用下继续慢慢下降，直到离月面 4 米高时再度悬停。此时，
关掉反冲发动机，探测器自由下落。由于探测器具备着陆缓冲机构，几个腿都有弹性，落地时不
至于摔坏。
安全降落以后，嫦娥三号将打开太阳能电池板接收能量，携带的仪器经过测试、调试后开始
工作。随后，“玉兔号”月球车将驶离着陆器，在月面进行 3 个月的科学勘测，着陆器则在着陆
地点进行原地探测。这将是中国航天器首次在地外天体的软着陆和巡视勘探，同时也是 1976 年
后人类探测器首次的落月探测。
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.1 题目要求
2. 嫦娥三号近月轨道示意图 (曲振东编制)
图 1.1: 嫦娥三号近月轨道示意图
®
3.
嫦娥三号着陆区域和着陆点示意图
图 1.2: 嫦娥三号着陆区域和着陆点示意图
¯
4. 主发动机和姿态调整发动机的分布图
嫦娥三号安装有大推力主减速发动机一台，位于正下方。小型姿态调整发动机 16 台，分布
在相对前、后、左、右四个侧面，如图 (1.3) 是一个侧面的分布情况。
图 1.3: 嫦娥三号主发减速动机与姿态调整发动机的分布图
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.1 题目要求 第一章 嫦娥 3 号
5. 关于比冲
比冲
°
或比冲量是对一个推进系统的燃烧效率的描述。比冲的定义为：火箭发动机单位质量
推进剂产生的冲量，或单位流量的推进剂产生的推力。比冲的单位为米/秒 (m/s)，并满足下列关
系式：
F
thrust
= v
e
˙m
其中：F
thrust
是发动机的推力，单位是牛顿；v
e
是以米/秒为单位的比冲； ˙m 是单位时间燃料消
耗的公斤数。
6. 关于月球参数
月球平均半径、赤道平均半径和极区半径分别为 1737.013km、1737.646km 和 1735.843km，
月球的形状扁率为 1/963.7256，月球质量是 7.3477 × 1022kg。月球与地球距离最远 (远地点)：
406610km，最近 (近地点)：356330km，平均距离为 384400km。
NASA 月球勘测轨道飞行器使用的月面海拔零点，是月球的平均半径所在的高度。所以，嫦
娥三号着陆点的海拔为 −2640m，即该点到月球中心的距离要比月球的平均半径少 2640m。
1.1.2 附件 2：嫦娥三号着陆过程的六个阶段及其状态要求；
1. 嫦娥三号软着陆过程示意图
图 1.4: 嫦娥三号软着陆过程示意图
°
维基百科，比冲. http://en.wikipedia.org/wiki/Specic_impulse. 2014.1.17
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.1 题目要求
2. 嫦娥三号软着陆过程分为 6 个阶段的要求
(1) 着陆准备轨道：着陆准备轨道的近月点是 15km，远月点是 100km。近月点在月心坐标
系的位置和软着陆轨道形态共同决定了着陆点的位置。
(2) 主减速段：主减速段的区间是距离月面 15km 到 3km。该阶段的主要是减速，实现到距
离月面 3 公里处嫦娥三号的速度降到 57m/s。
(3) 快速调整段：快速调整段的主要是调整探测器姿态，需要从距离月面 3km 到 2.4km 处
将水平速度减为 0m/s，即使主减速发动机的推力竖直向下，之后进入粗避障阶段。
(4) 粗避障段：粗避障段的范围是距离月面 2.4km 到 100m 区间，其主要是要求避开大的陨
石坑，实现在设计着陆点上方 100m 处悬停，并初步确定落月地点。
嫦娥三号在距离月面 2.4km 处对正下方月面 2300 × 2300m 的范围进行拍照，获得数字高程
如附图 (1.5) 所示 (相关数据文件见附件 3)，并嫦娥三号在月面的垂直投影位于预定着陆区域的
中心位置。
图 1.5: 距月面 2400m 处的数字高程图
该高程图的水平分辨率是 1m/像素，其数值的单位是 1m。例如数字高程图中第 1 行第 1 列
的数值是 102，则表示着陆区域最左上角的高程是 102 米。
(5) 精避障段：精细避障段的区间是距离月面 100m 到 30m。要求嫦娥三号悬停在距离月面
100m 处，对着陆点附近区域 100m 范围内拍摄图像，并获得三维数字高程图。分析三维数字高
程图，避开较大的陨石坑，确定最佳着陆地点，实现在着陆点上方 30m 处水平方向速度为 0m/s。
附图 (1.6) 是在距离月面 100m 处悬停拍摄到的数字高程图 (相关数据文件见附件 4)。
图 1.6: 距离月面 100m 处的数字高程图
该数字高程的水平分辨率为 0.1m/像素，高度数值的单位是 0.1m。
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.1 题目要求 第一章 嫦娥 3 号
(6) 缓速下降阶段：缓速下降阶段的区间是距离月面 30m 到 4m。该阶段的主要任务控制着
陆器在距离月面 4m 处的速度为 0m/s，即实现在距离月面 4m 处相对月面静止，之后关闭发动
机，使嫦娥三号自由落体到精确有落月点。
注：附件 3 和附件 4 中数字高程图对应的 *.tif 文件可以使用 Matlab 的“imread”命令打
开，“imread”的具体使用方法见 Matlab 相关帮助。
1.1.3 附件 3：距月面 2400m 处的数字高程图；
图 1.7: 距 2400m 处的数字高程图
1.1.4 附件 4：距月面 100m 处的数字高程图。
图 1.8: 距月面 100m 处的数字高程图
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
1.2 嫦娥 3 号软着陆的优化控制
1.2.1 基础知识准备
系统及最优控制基础知识
下面的内容介绍了一些系统和最优控制的基础知识。现代控制理论形成的主要标志是：1.Kalman
滤波理论；2.Pontryayin 极大值原理；3.Belman 动态规划方法。主要分支有线性系统理论、最优
控制理论、自适应控制、动态系统辨识、大系统理论等。
状态：是一些变量的集合，它们是足以完全描述系统全部运动的最小变量集。
完全描述：只要给定这组变量的某一初始时刻 t
0
的值，并且有输入控制函数 u(t)，那么系
统中的全部变量在 t ⩽ t
0
的值都唯一确定。
状态向量：如果
n
个状态变量用
x
1
(
t
)
, x
2
(
t
)
, . . . , x
n
(
t
)
表示，则
x(t) = [x
1
(t), x
2
(t), . . . , x
n
(t)]
T
记为状态变量。
状态方程：由系统的状态变量构成的一阶微分方程组称为状态方程。
示例：考虑下面的例子，线性弹簧-质量-阻尼器构成的机械位移系统 (k − m − f ) 如图 (1.9)
所示
图 1.9: 弹簧阻尼位移系统
如果输入量是外力 F (t)，输出量是位移 y(t)，则
m
d
2
y
dt
2
+ f
dy
dt
+ ky = F (t)
y(t),
dy
dt
就是足以完全描述系统运动的最小一组变量。输入量 u(t) = (F
t
, θ
t
)，输出量 y(t) = m
t
，
状态量 x(t) = (v, r, w)。
上面介绍了系统的一些基本定义，下面我们介绍基于系统的最优控制问题，以及用变分法、
P 极大值原理、动态规划等技术求解系统的最优控制问题。常见的最优控制包括最短时、最少能
耗及其组合优化问题。
经典控制论中常用一个高阶微分方程来描述系统的运动规律，而现代控制论中采用是状态
空间法，即为一组状态变量的一阶微分微分方程的形式作为系统的数学模型。这使得我们可以求
解系统中许多量的值。MATLAB：控制系统工具箱 Control System Toolbox、系统辨识工具箱
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
System Idemtication Toolbox、鲁棒性控制工具箱 Robust Control Toolbox、多变量频域设计工
具箱 Muli-Vanable Frequency Design Toolbox、最优化工具箱 Ootimization Control toolbox 以
及 Similink 控制系统仿真环境。
抽象系统方程为
˙x(t) = F (x(t), u(t), t)
其中：x(t) 为 n 个状态变量 x
1
, x
2
, . . . , x
n
的向量形式；u(t) 为控制输入 u
1
, u
2
, . . . , u
r
的向量形
式；F = (f
1
, f
2
, . . . , f
n
)
T
是一个函数标量。
输出方程为
Y (t) = G(x(t), u(t), t)
其中：G = (g
1
, g
2
, . . . , g
m
)
T
是函数矢量。
现代控制理论中，用系统状态方程和输出方程来描述系统的状态行为，状态方程和输出方程，
合起来称为系统的状态表达式/动态方程。根据一阶微分方程 (F, G) 的形式，一般系统可分为：1.
线性定常系统；2. 线性不定常系统；3. 非线性定常系统；4. 非线性不定常系统。即
˙x = Ax
˙x = A(t)x
˙x = Ax
2
+ x
˙x = A(t)x
2
+ x
对于一个系统的分析，一般有以下步骤：1. 系统动态方程表达；2. 系统状态方程的求解；3. 系统
稳定性及李雅普诺夫稳定；4. 能控制和能观性；5. 系统的最优控制。
在最优控制过程中，我们不仅要给出系统的动态方程，而且还要注意系统变量的额约束。对
系统状态 x(t) 而言，要考虑系统的始、终端条件。始终端约束条件即系统在时间 t
0
和 t
f
时，状
态 x, x
f
应该满足一定的条件。端点条件一般有 3 种类型：1. 固定端；2. 自由端；3. 可变端 (灰
色端)。
固定端就是时间 t
0
, t
f
固定，状态值 x
0
, x
f
固定；自由端是指时间 t
0
, t
f
固定，但端点状态
x
0
, x
f
不变受任何限制的端点；可变端是指 t
0
, t
f
, x
0
, x
f
都变的端点，但其一般应该满足某些条
件，例如 N[x
f
, t
f
] = 0，灰色端 ∆
1
⩽ x
f
⩽ ∆
2
。
对系统控制变量 u
t
而言，控制输入 u(t) 往往不能不受限制，因为其往往表示非常实际的物
理量 (例如：推理 F
t
，夹角 θ
t
)。因此，实际中其约束一般表示 u
t
在某一范围内取值，如下表示
u
min
⩽ u(t) ⩽ u
max
由控制约束条件所规定的点集称为控制域，记为 R
u
。
凡在区间 [t
0
, t
f
] 内定义，且在 R
u
内取值的每一个控制函数 u
t
均称为容许控制，u(t) ∈ R
u
。
最优控制的目标泛函：系统的性能指标一般都是一个控制信号 ( ˙m(t) 函数) 的函数，即泛函。
对连续时间系统，目标泛函一般为
J = Φ[x(t
f
)] +
t
f
t
0
L(x(t), u(t), t)dt
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
上述目标形式称为综合型/Bolza 型。
J = Φ[x(t
f
)]
称为终端型/Mayer 型。
J =
t
f
t
0
L[x(t), u(t), t]dt
称为积分型/Langrange 型。其中：L 为标量函数。最优控制问题即是在问题空间中，求解最优控
制 u(t)
∗
，使目标泛函 J 最优。
解法：最优控制问题实际上就是微分代数混合优化问题，一般有两种求解方法：解析法和数
值法。解析法又分为变分法和贝尔曼最佳原理的动态规划方法，但我们所建立的许多模型都相对
复杂，解析法并不是非常适用，所以在书中大多数地方，我们都讨论问题的数值解法。下面，我
们将简单介绍一下经典的最优控制解析方法：1. 变分法；2.P 极大值原理；3. 动态规划。
变分法
考虑无约束情况下的极值问题。(1) 一元函数极值：设 f(x) 为定义在 [a, b] 上的连续可微函
数，且二阶可微。如果
dy
dx
x=x
0
= f
′
(x)|
x=x
0
= 0
d
2
y
dx
2
x=x
0
= f
′′
(x)|
x=x
0
> 0
则 f 在 x
0
处为极小值。
(2) 二元函数极值：设 f(x, y) 为定义在 D 上的连续可微函数，且二阶可微。如果
∂f
∂x
x=x
0
∂f
∂y
y=y
0
且
∂
2
f
∂x
2
< 0，则 f 在点 (x
0
, y
0
) 处获得极小值。
(3) 多元极值：对一个 n 元函数 J = f (x)，取极值的必要条件
∂f
∂x
= 0
取极小值的充分必要条件是
∂
2
f
∂x
2
> 0
其中：
∂
2
f
∂x
2
= H 为 Hessi 矩阵，若 Hessi 矩阵为负定，则函数存在极大值。
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
(4) 泛函的定义：泛函可以简单理解为函数的函数。设对自变量 t 存在一类函数 {x(t)}，如
果 ∀x
t
∈ ϕ，有一个 J 值与之对应，则变量 J 称为依赖于函数 x(t) 的泛函数，简称泛函，如图
(1.10) 所示
图 1.10: 泛函示意图
将 {x
i
(t)}
4
i=1
带入某关系式，有 J = Af [x(t)]。故称 J 为依赖 x(t) 的泛函，它是普通函数
的一种扩充，函数 x(t) 称为 J 的宗基。当然，{x(t)} 可以是一参数函数族：x(t|θ)。
由上述泛函定义可见，泛函是一标量，例如
J[x(t)] =
1
0
x(t)dt
当 x(t) =
1
2
t 时，J = 1；当 x(t) = cos t 时，J = sin 1。在最优控制中
J =
t
f
t
0
L[x(t), u(t), t]dt + ∆
J
取决于
u
(
t
)
, x
(
t
)
，故
J
为泛函。
在学习函数时，我们讨论了函数的连续性、函数导数和极值问题，那么，对于泛函，我们仍
然讨论与之相似的问题。
(5) 泛函的连续性：称宗量 x(t) 的增量 (两函数之差) 为 x(t) 的变分
δx = x(t) − x
0
(t)
那么有
δ ˙x =
d
dx
(δx)
即导数的变分等价于变分的导数，且 δx 是变量 t 的函数 δx(t)。
定义 (泛函的连续性) 如果对于 x(t) 的微小变化，有 J[x(t)] 的微小变化与之对应，则说
J[x(t)] 是连续的。对 ∀ϵ > 0，∃δ > 0，当
|x(t) − x
0
(t)| < δ
| ˙x(t) − ˙x
0
(t)| < δ
.
.
.
| ˙x
k
(t) − x
k
0
(t)| < δ
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
时，有 |J[x(t)] − J[x
0
(t)]| < ϵ，则称 J[x(t)] 在 x
0
(t) 处是 k 阶接近连续的。
泛函的连续性如图 (1.11) 所示
图 1.11: 泛函的连续性示意图
(6) 泛函的变分：在函数中，微分是函数增量 ∆y 对自变量 ∆x 的线性部分。那么，在泛函
中亦可将变分视为泛函增量 ∆J 对宗量 x(t) 的线性部分，即当宗量 x(t) 有增量 δx(t) 时，泛函
J[x(t)] 的增量可以表示为
∆J = J[x(t) − δx(t)] − J[x(t)] = L[x(t), δx(t)] + r[x(t), δx(t)]
其中：L[x(t), δx(t)] 是泛函增量的线性主部分。故
δJ = L[x(t), δx(t)]
称为 J 的变分。
泛函变分的求法：J[x(t)] 的变分等于泛型 J[x(t) + αδx(t)] 对 α 的导数在 α = 0 时的值
δJ =
∂
∂α
J[x(t) + αδx(t)]
α=0
= L[x(t), δx(t)]
(7) 泛函极值：泛函极值的定义
定义 (泛函极值) 如果 J 在任意一点与 x
0
(t) 接近的曲线 x(t) 上的值不小于 J[x
0
(t)]，即
J[x(t)] − J[x
0
(t)] ⩾ 0
则 J 在 x
0
(t) 上达到极小值。
泛函极值的必要条件 (无约束)：若 J 在 x
0
(t) 处达到极小值，则在 x(t) = x
0
(t) 上有
δJ = 0
上述欧拉方程是泛函极值的必要条件而非充分。
固定端点的泛函极值问题：已知宗量 x(t) 在 t
0
时为 x
0
，在 t
f
时为 x
f
，则使积分型性能指
标函数
J =
t
f
t
0
L[x(t), ˙x(t), t]dt
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
取极值的必要条件为 x
∗
(t) 满足下面的欧拉方程
∂L
∂x
−
d
∂t
∂L
∂ ˙x
= 0
其中：L 及 x(t) 在 [t
0
, t
f
] 上至少两次连续可微。
将上述结论推广至多变量系统
J
=
t
f
t
0
L[x(t), ˙x(t), t]dt
x(t
0
) = x
0
x(t
f
) = x
f
J 极值存在的必要条件为
∂L
∂x
−
d
∂t
∂L
∂ ˙x
= 0
及横截条件方程
∂L
∂ ˙x
T
t=t
f
δx(t
f
) −
∂L
∂ ˙x
T
t=t
0
δx(t
0
) = 0
其中：
∂L
∂x
=
∂L
∂x
1
,
∂L
∂x
2
, . . . ,
∂L
∂x
n
T
∂L
∂ ˙x
=
∂L
∂ ˙x
1
,
∂L
∂ ˙x
2
, . . . ,
∂L
∂ ˙x
n
T
示例：求泛函
J(x
1
, x
2
) =
π
2
0
[2x, x
2
+ ˙x
2
1
+ ˙x
2
2
]dt
x
1
(0) = 0, x
2
(π/2) = 1
x
2
(0) = 0, x
2
(π/2) = 1
下的极值。求解程序如下
1 [ x1 , x2 ] = ds ol ve ( ’D2x1−x2=0 ’ , ’D2x2−x1=0 ’ , . . .
2 ’ x1 ( 0)=0 ’ , ’ x1 ( pi /2)=1 ’ , ’ x2 ( 0)=0 ’ , ’ x2 ( p i /2)=−1’ , ’ t ’ )
3
可变端点极值问题的必要条件：设宗量 x(t) 自一给定点 (t
0
, x
0
) 达到给定 x(t
f
) = φ(t
f
) 上
某一点 [t
f
, φ(t
f
)]，使
J =
t
f
t
0
L[x(t), ˙x(t), t]dt
取极值的必要条件为 x
∗
(t) 满足如下欧拉方程
∂L
∂x
−
d
∂t
∂L
∂ ˙x
= 0
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
及
x(t
0
) = x
0
L[x, ˙x, t] + ( ˙φ − ˙x)
T
∂L
∂ ˙x
t=t
f
= 0
x(t
f
) = φ(t
f
)
其中：x(t) 连续两阶可导；L 至少两次连续可微；φ(t) 有连续一次导数。可变端点极值问题示意
图如图 (1.12) 所示
图 1.12: 可变端点极值问题示意图
变分法求解最优控制问题
前面导出了极值曲线 x
∗
(t) 存在的必要条件-欧拉方程和横截条件 (x(t) 无约束)，而在最优
控制中，x(t) 除了要满足限制条件外，还应该满足某些系统动态方程，它可以看成一种等式约束，
可采用拉格朗日乘子法进行处理。将等式约束下对 J 的优化转化为无约束条件下就 Hamilton 的
极值问题。它只适用于控制变量无约束的情况。
min
u(t)
J =
t
f
t
0
L[x(t), u(t), t]dt
s.t.
˙x(t) = f [x(t), u(t), t]
x(t
0
) = x
0
x(t
f
) = x
f
¬固定端点的控制问题；可变端点的控制(t
f
固定，x(t
f
) 自由；t
f
固定，x(t
f
) 有约束；t
f
可
变，x(t
f
) 有约束)。对问题¬(1.2.1)，引入拉格朗日乘子 λ(t) = [λ
1
(t), λ
2
(t), . . . , λ
n
(t)]
T
，构建哈
密顿标量函数
H[x(t), u(t), λ(t), t] = L[x(t), u(t), t] + λ
T
(t)f[x(t), u(t), t]
则上述优化问题变为
J =
t
f
t
0
H[x(t), u(t), λ(t), t] − λ
T
x(t)dt
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
由欧拉方程得
∂[H − λ
T
˙x]
∂x
−
d
dt
∂[H − λ
T
˙x]
∂ ˙x
= 0
⇒
˙
λ = −
∂H
∂x
伴随方程
∂[H − λ
T
˙x]
∂u
−
d
dt
∂[H − λ
T
˙x]
∂ ˙u
= 0
⇒
∂H
∂u
= 0 控制方程
∂[H − λ
T
˙x]
∂λ
−
d
dt
∂[H − λ
T
˙x]
∂
˙
λ
= 0
⇒ ˙x = −
∂H
∂λ
= f(x, u, t) 状态方程
及边界条件 (横截条件)
∂[H − λ
T
˙x]
∂ ˙x
t=t
f
= 0
⇒ λ(t
f
) = 0
状态方程与伴随方程合称为正则方程。
求解过程：
∂H
∂u
= 0 可求得最优控制 u
∗
与 x 及 λ 的函数关系式，然后将其带入正则方程组
中消去 u
∗
，解两点边值问题
˙x = f (x, u(x, λ), t), x(t
0
) = x
0
˙
λ = −
∂H
∂x
, x(t
f
) = x
f
可求得 x
∗
(t), λ
∗
(t)，再将 x
∗
, λ
∗
带入 u(x
∗
, λ
∗
) 可得 u
∗
的具体情况。另外，对 H 函数有
dH
d
t
=
∂H
∂t
+
∂H
∂u
T
˙u +
∂H
∂x
T
˙x +
∂H
∂λ
T
˙
λ
=
∂H
∂t
+
∂H
∂u
T
˙u +
∂H
∂x
+
˙
λ
T
f
当
∂H
∂u
∂H
∂x
= λ,
∂u
∂λ
= 0 时，
dH
dt
=
∂H
∂t
，即沿最优曲线，H 对 t 的导数等于对时间的偏导数。注：
上述欧拉方程式取极值的必要条件。
下面，分析终端自由的情况(1.2.1)：
J = Φ[x(t
f
), t
f
] +
t
f
t
0
L(x, u, t)dt
s.t.
˙x = f (x, u, t)
t
0
, x
0
固定
求 t
f
, x
f
自由时的最优控制 u
∗
及最优轨线 x
∗
。
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
(1) 构建 H 函数
H(x, u, λ, t) = L(x, u, t) + λ
T
f(x, u, t)
(2) 增广泛函 J
′
J
′
= Φ[x(t
f
), t
f
] +
t
f
t
0
[H − λ
T
x]dt
J
′
的一次变分为
δJ
′
=
∂Φ
∂x(t
f
)
T
δx(t
f
) +
∂Φ
∂x(t
f
)
T
∂x(t
f
)
∂t
f
δt
f
+
∂Φ
∂t
f
δt
f
+ [H − λ
T
˙x]
t
f
δt
f
+ δJ
0
其中：δJ
0
= δJ
0
[x, ˙x, λ, u] = δ
t
f
t
0
[H − λ
T
˙x]dt。
泛函 J
′
极值存在的必要条件是：δJ
′
= 0。对固定时间 t
0
, t
f
的变分，δ(t
0
) = 0
δJ
0
=
t
f
t
0
∂
¯
H
∂x
−
d
dt
∂
¯
H
∂ ˙x
T
δxdt +
t
f
t
0
∂
¯
H
∂u
T
δudt
+
t
f
t
0
∂
¯
H
∂λ
T
δλdt +
∂
¯
H
∂
˙
x
T
δx
t
f
t
0
其中：
¯
H = H − λ
T
x。将 δJ
0
带入 δJ
′
，且令 δJ
′
= 0，有
˙
λ = −
∂H
∂x
˙
x
=
−
∂H
∂λ
正则方程
∂H
∂u
= 0 控制方程
边界条件/横截条件
λ(t
f
) =
∂Φ
∂x
t=t
f
x(t
0
) = x
0
H[x(t
f
), u(t
f
), t
f
] = −
∂Φ
∂t
f
¬终端 t
f
固定，x(t
f
) 自由时
x(t
f
) =
∂Φ
∂x
t=t
f
x(t
0
) = x
0
终端 t
f
固定，x(t
f
) 有约束：N
1
(x(t
f
), t
f
) = 0，其中：N
1
= [N
11
, N
12
, . . . , N
1m
]
T
λ(t
f
) =
∂Φ
∂x
+
∂N
T
1
∂x
v
t=t
f
x(t
0
) = x
0
N
1
(x(t
f
), t
f
) = 0
http://www.ma-xy.com 15 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
®t
f
可变，x(t
f
) 有约束：N
1
(x(t
f
), t
f
) = 0
λ(t
f
) =
∂Φ
∂x
+
∂N
T
1
∂x
v
t=t
f
H +
∂Φ
∂t
+ v
T
∂N
1
∂t
t=t
f
= 0
x(t
0
) = x
0
N
1
(x(t
f
), t
f
) = 0
¯当 t
f
, x(t
f
) 稳定时
dx(t
f
) = x(t
f
+ δt
f
) + δx(t
f
+ δt
f
) − x(t
f
)
δx(t
f
) + ˙x(t
f
)δt
f
°当 x
f
固定，t
f
自由时
dx(t
f
) = 0 = δx(t
f
) + ˙xδt
f
极大 (小) 值原理求解最优控制问题
上面介绍的变分法用于求解 u(t) 无约束的情况，而在实际过程中 u(t) 往往具有实际的物理
意义，因此 u(t) 一般有约束，例如
A ⩽ u(t) ⩽ B
设 R = {u(t)|∀t ∈ [t
0
, t
f
], A ⩽ u(t) ⩽ B} 为函数空间。P 极大值原理用于求解 u(t) 有约束的情
况
min
u(t)
J = Φ[x(t
f
), t
f
] +
t
f
t
0
L[x(t), u(t), t]dt
s.t.
x(t
0
) = x
0
˙x = f (x(t), u(t), t)
N
1
[x(t
f
), t
f
] = 0
g[x(t), u(t), t] ⩾ 0
如何处理上面的 g ⩾ 0？要是能将 g 变为等式约束就好了，我们仿照线性规划，引入人工变
量 z, w
( ˙z)
2
= g[x, u, t], z(t
0
) = 0
˙w = u(t), w(t
0
) = 0
无论 ˙z 是正还是负，( ˙z)
2
恒非负。如 u(t) 是分段连续函数免责 w(t) 是分段光滑连续函数。如此，
将 g ⩽ 0 转化为 ( ˙z)
2
= g[x(t), ˙w(t), t] 等式约束，继而可以引入拉格朗日乘子，构建 H 函数进行
求解。
http://www.ma-xy.com 16 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
实现最优控制的必要条件是：
(1)u
∗
, x
∗
, λ
∗
满足下面条件
˙x =
∂H
∂λ
˙
λ = −
∂H
∂x
其中：H = L + λ
T
f。
(2)
H[x
∗
, u
∗
, λ
∗
, t] = min
u(t)∈R
H[x
∗
, u, λ
∗
, t]
(3) 端点边界条件与终端横截条件
x(t
0
) = x
0
N[x(t
f
), t
f
] = 0
λ
f
=
∂Φ
∂x
+
∂H
1
∂x
T
v
t=t
f
(4) 终端 t
f
可变时
H +
∂Φ
∂t
+ v
T
∂N
1
∂t
t=t
f
= 0
最小时间控制
J = t
f
砰-砰
最小能量控制
J =
t
f
t
0
|u(t)|dt 砰-零-砰
最小时间能量综合控制
J =
t
f
t
0
ρ + |u(t)|dt 砰-零-砰
其中：ρ 为权重。
最优控制问题的数学提法
用极大值原理结合非线性规划来解决软着陆最小消耗问题，具有算法简单、原理可靠的特点，
但其只能使发动机推力 F 为常量，且由于末端时刻未定，导致迭代次数较大。此外，还需要猜测
无物理意义的协变量状态的初始值，这使系统很敏感。用“协状态-控制”变换可以从一定程度上
减弱这种对初值的依赖性。为了避免上述缺点，下面引入参数化控制的方法，该方法是 K.L.Teo
和 H.E.Jlee 等人提出，主要思想是用若干个分段常值函数逼近最优解。
http://www.ma-xy.com 17 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
(1) 设控制系统的运动遵从抽象微分方程
˙x = f (x(t), u(t), t) t ∈ J
其中：f ∈ [X × U × J; Y ]；J 为 [0, ∞) 中的某个连续闭子集，称为系统的时间域；X 为 Banach
空间/拓扑线性空间，称为系统的相空间或状态空间；U 为拓扑空间，称为系统的控制域；t 为时
间变量；x 为状态变量；u 为控制变量；x 的具体情况为状态轨线，是时间 t 的函数，其图形是
J × X 中的点集；相空间 X 中确定的点集 {x(t)|t ∈ J} 称为系统的一条相轨线，与 t 无关。
(2) 控制函数 u 的限制。例如要求 u(·) 为简单函数、可测函数、连续函数、分段函数、分段
可微等，即
{u} = U ⊂ [J; U]
(3) 端点条件。有 t
0
=
J, t
f
= sup J，对于某个给定的点集 M = M
0
× M
1
⊂ Y × X
x(t
0
) ∈ M
0
x(t
f
) ∈ M
1
当 M
1
= X 时，称为终端自由；当 t
f
= sup J < ∞ 时，称终端时间 tf 固定。在许多控制问题
中，控制系统的状态空间 X 是某个广义函数空间或抽象的向量值函数空间。
假设一个控制系统的状态方程、控制域 U 、控制函数类 U 以及约束条件均已给出，且相应
的容许控制类 U
ad
= 0。记全体容许对所构成的集合为 A
ad
，任一映射 J ∈ [A
sd
; R] 称为该控制
系统的一个指标泛函。
对于给定的指标泛函 J，寻找适当的 (x, u)
∗
∈ A
ad
，使得
J(x, u)
∗
= J
∗
= inf{J(x, u)|(x, u)
∗
∈ A
ad
}
如果这样的容许对 (x, u)
∗
存在，则称该最优控制问题有解。(x, u)
∗
为最有对，u
∗
为最优控制，
x
∗
为最优轨线。
给出非线性标量型的最优化问题
J = ϕ
0
[x(t
f
|u)] +
t
f
t
0
L
0
[x(t|u), u(t), t]dt
s.t. g
i
(u) = ϕ
i
[x(t
f
|u)] +
t
f
t
0
L
i
[x(t|u), u(t), t]dt ⩾ 0
其中：Q
i
, L
i
(i = 1, 2, . . . , N ) 实值向量函数；f : [0, t
f
] × R
n
s
× R
n
c
→ R
n
s
；ϕ
i
: R
n
s
→ R, i =
0, 1, . . . , N,L
i
= [0, t
f
] × R
n
c
× R
n
c
→ R, i = 0, 1, . . . , N。
∃k，∀(t, x, u) ∈ [0, t
f
] × R
N
S
× V ，均有 |f (t, u, x)| ⩽ k(1 + |x|)。其中：V ∈ R
n
c
是 R
n
c
中
任意紧致闭子集。泛函 f 和 L
i
, i = 0, 1, . . . , N ，关于 x 和 u 的偏导数在 [t
0
, t
f
] 内均为分段连续
函数。∀i, ϕ
i
对 x, u 是连续可微的。
参数化最优控制器设计。首先，选取一组满足 {t
p
k
}
n
p
k=0
t
0
= t
p
0
< t
p
1
< t
p
2
< · · · < t
p
n
p
= t
f
http://www.ma-xy.com 18 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
和一组参数 σ
p
k
, k
= 1
,
2
, . . . , n
p
，然后构造形如
u
p
(t) =
n
p
i=1
σ
p
k
χ
[t
p
k−1
,t
p
k
)
(t)
所示的参数化分段常数控制器。其中：
χ
[t
p
k−1
,t
p
k
)
(t) =
1 t ∈ [t
p
k−1
, t
p
k
)
0
令 σ
p
= [σ
p
1
, . . . , σ
p
n
p
]，将控制器带入系统方程可得
˙x(t) =
˜
f[x(t), t, σ
p
]
= f
x(t), t,
n
p
i=1
σ
p
k
χ
[
t
p
k
−
1
,t
p
k
)
(t)
于是，最优控制问题变为寻找一组最优参数 σ
p
，最小化指标泛函
J(σ
p
) = ϕ
0
[x(t
f
|σ
p
)] +
t
f
t
0
˜
L
0
[x(t, σ
p
), σ
p
, t]dt
s.t.
g
i
(
σ
p
) =
ϕ
i
[
x
(
t
f
|
σ
p
)] +
t
f
t
0
˜
L
i
[x(t|σ
p
), σ
p
, t]dt ⩽ 0
˜
L
i
[x, σ
p
, t] = L
i
[· · · ]
i = 0, 1, 2, . . . , N
显然，对于每个给定的 p，这都是一个有限维的参数化问题，且当 p → ∞ 时，最优解收敛。
指标泛函和约束条件关于参数 σ
p
的梯度公式为
∂J
∂σ
p
=
t
f
t
0
∂
˜
H
0
[t, x, σ
p
, λ
0
(σ
p
)]
∂σ
p
dt
∂g
i
∂σ
p
=
t
f
t
0
∂
˜
H
i
[t, x, σ
p
, λ
i
(σ
p
)]
∂σ
p
dt
其中：
˜
H
i
[t, x, σ
p
, λ
i
(t|σ
p
)] =
˜
L
i
[t, x, σ
p
, λ
i
(t|σ
p
)]+ λ
T
˜
f[t, x(t|σ
p
), σ
p
]；λ = [λ
1
, . . . , λ
n
s
]
T
；
˙
λ
i
(t|σ
p
) =
−
∂
˜
H
∂x(t|σ
p
)
。
边值条件：
λ
i
(t
f
) =
∂
˜
ϕ
i
[x(t
f
|σ
p
)]
∂x(t
f
|σ
p
)
上面所提到的参数化控制器只将 σ
p
视为参数，t
k
取定值，而 t
k
取值不同会影响到参数控
制器的逼近程度。相关研究表明，在数值计算中，如果直接将其视为参数，则求解参数梯度时难
度很大，甚至不能求解，为此我们引入强化技术。
强化技术：从 s ∈ [0, 1] 到 t ∈ [t
0
, t
f
]，构建如下变换
t(s) =
δ
p
k
(s − ξ
p
k−1
) k = 1
k−1
j=1
δ
p
j
(ξ
p
j
− ξ
p
j−1
) + δ
p
k
(s − ξ
p
k−1
) k = 2, 3, . . . , n
p
http://www.ma-xy.com 19 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
其中：δ
p
k
>
0
，
s
∈
[
ξ
p
k−1
, ξ
p
k
]，{ξ
p
k
}
n
p
k=1
为 [0, 1] 区间上预先给定的分段点，且满足 0 = ξ
p
0
< · · · <
ξ
p
k−1
< ξ
p
k
< · · · < ξ
p
n
p
= 1。
将上式两边对 s 求导，对不可导的点用左导数替代常规导数，有
dt(s)/ds = v
p
(s)
其中：v
p
(s) =
n
p
k=1
δ
p
k
χ
[ξ
p
k−1
,ξ
p
k
]
(s)；当 s ∈ [ξ
p
k−1
, ξ
p
k
) 时，χ
[ξ
p
k−1
,ξ
p
k
)
(s) = 1，否则为 0，即 χ 为特征
函数。不妨令
ˆx(s) =
x(t(s))
t(s)
ˆu(s) =
u
p
(t(s))
v
p
(s)
ξ
p
= [ξ
p
1
, . . . , ξ
p
n
p
]
将其带入到 ˙x(t) =
˜
f(t, x, σ
p
)，得到增广系统
dˆx(s)
ds
=
v
p
(s)
ˆ
f(s, ˆx(s) u
p
(s))
v
p
(s)
其中：
ˆ
f(s, ˆx(s), u
p
(s)) = f (t(s), x(t(s)), u
p
(t(s)))。系统的初始条件为：˜x(0) = x
0
, t(0) = 0。终
端条件为 t(1) = t
f
，则最优控制为
J(σ
p
, r
p
) = ϕ
0
(˜x(1|σ
p
, r
p
)) +
1
0
ˆ
L
0
(t(x), ˜x(s|σ
p
, r
p
), σ
p
, r
p
)ds
s.t.
g
i
(σ
p
, r
p
) = ϕ
i
(˜x(1|σ
p
, r
p
)) +
1
0
L
i
(t(s), ˜x(s|σ
p
, r
p
), σ
p
, r
p
)ds ⩾ 0
ˆ
L
i
(t(x), ˜x(s|σ
p
, r
p
), σ
p
, r
p
) = v
p
(s|r
p
)
˜
L
i
(t(s), x(t(s)|σ
p
), σ
p
)
及边值条件。其中：ξ
p
k
为时间分段点。指标泛函及约束函数关于 r
p
求导有
∂J
∂r
p
=
∂ϕ
0
(˜x(1|σ
p
, r
p
))
∂r
p
+
1
0
∂
ˆ
H
0
(t(s), ˜x(s|σ
p
, r
p
), σ
p
, r
p
, λ
0
(s|σ
p
, r
p
))
∂r
p
ds
∂g
i
∂r
p
=
∂ϕ
i
∂r
p
+
1
0
∂H
i
(−λ
i
)
∂r
p
ds
其中：λ
i
= [λ
it
, . . . , λ
ix
]
ˆ
H
i
(t, ˜x, σ, r, λ
i
) =
ˆ
L
i
+ λ
nt
v
p
+ λ
ix
v
p
f(t, ˜x, σ, r)
协状态变量 λ
i
(s|σ
p
, r
p
), s ∈ [0, 1] 可由下列方程解的
d(λ
ix
(s))
ds
= −
∂
ˆ
H
i
(t, ˜x, σ
p
, r
p
, λ
i
)
∂ ˜x
dλ
it
(s)
ds
= −
∂
ˆ
H
i
∂t(s)
http://www.ma-xy.com 20 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
边值条件
λ
ix
(1) =
∂ϕ
i
(˜x(1|σ
p
, r
p
))
∂ ˜x
λ
it
(1) =
∂ϕ
i
(˜x(1|σ
p
, r
p
))
∂t
约束变换技术：h(t, x, u) ⩾ 0 的约束条件是常见的，但我们要将其转化为规范标准形式，当
h(t, x, u) = 0 时
g(u) =
t
f
0
h(t, x, u)
2
dt = 0
当 h ⩾ 0 时，Teo 于 1993 年给出了一种约束变换技术
g(u) =
t
f
0
min{h(t, x, u), 0}dt = 0
然而，上式在 h = 0 处并不光滑，因此用下面方法去接近
g =
t
f
0
L
ε
(h)dt + τ ⩾ 0
其中：
L
ε
(h) =
h h ⩾ ε
(h − ε)
2
/4ε −ε < h < ε
0 h ⩽ −ε
ε, τ > 0 是调节参数，当 ε 足够小或 τ(ε) > 0，使得 ∀τ, 0 < τ < τ (ε)，能够令不等式达到满足要
求。
对原始问题引用约束变换技术，参数化控制器及强化技术，有
dˆr(s)
ds
= v
p
(s)ˆv(s)
dˆv(s)
ds
= v
p
(s)
ˆ
F (s)
ˆm(s)
sin
ˆ
ψ(s) −
G
0
M
ˆr
2
(s)
+ ˆr(s) ˆw
2
(s)
d
ˆ
θ(s)
ds
= v
p
(s) ˆw(s)
d ˆw(s)
ds
= v
p
(s)
−
ˆ
F (s)
ˆm(s)ˆr(s)
cos
ˆ
ψ(s) +
2ˆv(s) ˆw(s)
ˆr(s)
d ˆm(s)
ds
= −v
p
(s)
ˆ
F (s)
v
e
dt(s)
ds
= v
p
(s)
http://www.ma-xy.com 21 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
初始条件为：ˆr(0) = r
0
, ˆv(0) = 0,
ˆ
θ(0) = 0, ˆw(0) = w
0
, ˆm(0) = m
0
, t(0) = 0。和控制器
ˆu(s) =
ˆ
F
p
(s) =
n
p
i=1
σ
p
F,i
χ
[ξ
p
i−1
,ξ
p
i
)
(s)
ˆ
ψ
p
(s) =
n
p
i=1
σ
p
ψ,i
χ
[ξ
p
i−1
,ξ
p
i
)
(s)
ˆv
p
(s) =
n
p
i=1
δ
p
i
χ
[ξ
p
i−1
,ξ
p
i
)
(s)
取参数 {σ
p
F,i
}
n
p
i=1
, {σ
p
ψ,i
}
n
p
i=1
, {δ
p
i
}
n
p
i=1
，使
min
ˆ
J = ˆm(0) − ˆm(1)
s.t.
ˆr(1) = r
f
ˆv(1) = v
f
ˆw(1) = 0
g =
t
f
0
v
p
Lε(ˆr)dt + τ > 0
L
ε
(ˆr)
求解算法：
Step1. 给定分点个数 n
p
，约束变换技术超参 ε, τ ，选定区间 [0, 1] 区间的单调序列 {ξ
p
k
}
n
p
k=1
。
Step2. 任选一组参数 {σ
p
F,i
}
n
p
i=1
, {σ
p
ψ,i
}
n
p
i=1
, {σ
p
i
}
n
p
k=1
Step3. 将参数带入系统，并求取指标函数关于参数的梯度。
Step4. 利用经典的优化算法更新 {σ
p
F,i
}, {σ
p
ψ,i
}
,
{σ
p
i
}。
Step5. 判断更新后的参数是否满足终端条件和不等式约束，同时，
ˆ
J 的值是否满足要求，不满
足则返回 Step2.
1.2.2 模型假设
1. 假设月球是一个球体；
2. 假设只有月球对飞行器有影响；
3. 忽略月球自转；
4. 嫦娥三号软着陆轨道为过月球自转轴的平面；
1.2.3 符号说明
1.2.4 问题一的分析与求解
问题的分析
确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。
http://www.ma-xy.com 22 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
模型的建立与求解
为方便，我们假设月球是一个球体。飞行器 (嫦娥 3 号) 绕月飞行如图 (1.13) 所示
图 1.13: 绕月飞行轨道示意图
其中：v
A
为近月点的速度，v
B
为远月点速度，r 为月球半径，d
1
为近月点到月面的距离，d
2
为
远月点到月面的距离。
假设只有月球对飞行器有影响。下面，我们先来求近月点和远月点的速度 v
A
, v
B
。
(1) 机械能守恒定律。由于在绕月飞行阶段，飞行器不消耗燃料，故有机械能守恒：机械能
= 动能 + 势能。势能可以是弹性势能/重力势能。只在重力做功的情况下，物体的动能和势能之
间相互转换，但总和不变
E
A
=
1
2
mv
2
A
−
GMm
r
A
E
B
=
1
2
mv
2
B
−
GMm
r
B
其中：E 为机械能，M 为月球质量，m 为飞行器质量，G 为万有引力常量，r
A
= r + d
1
为近月
点到月心的。由机械能守恒，我们有
E
A
= E
B
(1.1)
(2) 开普勒第二定律。1609 年《新天文学》德国天文学家约翰尼斯. 开普勒指出：在相等时
间内，太阳和绕太阳运行的行星的连线所扫过地面积是相等的。并且该定律适用于宇宙中一切绕
心的天体运动。单位时间内飞行器在近月点和远月点扫过的面积为
S
A
=
1
2
r
A
v
A
S
B
=
1
2
r
B
v
B
由开普勒第二定律，我们有
S
A
= S
B
(1.2)
http://www.ma-xy.com 23 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
我们要求两个参数 v
A
, v
B
，就需要两个方程，组合方程 (1.1) 和方程 (1.2) 有
E
A
= E
B
S
A
= S
B
⇒
v
A
=
2GMr
B
r
A
(r
A
+ r
B
)
= 1692.7m/s
v
B
=
2GMr
A
r
B
(r
A
+ r
B
)
= 1614.4m/s
其中：M = 7350 × 10
22
kg，G = 6.672 × 10
−11
Nm
2
/kg
2
，r
A
= r + d
1
= 1752013m，r
B
= r + d
2
=
1837013m。
(3) 开普勒第一定律。除了上面的这种方法之外，我们还可以寻找其他方程来求解速度 v
A
, v
B
。
下面，我们来介绍开普勒第一定律。开普勒在《宇宙和谐论》中表述：每一个行星都沿各自的椭
圆轨道环绕太阳，而太阳则处于椭圆的一个焦点上。如图 (1.14) 所示
图 1.14: 开普勒第一定律示意图
在图 (1.14) 中的近月点 A 和远月点 B 处，我们有
GMm
r
2
A
= m
v
2
A
r
r =
b
2
a
b
2
= a
2
− c
2
2a = ar
0
+ d
1
+ d
2
GMm
r
2
B
= m
v
2
B
r
r =
b
2
a
b
2
= a
2
− c
2
2a = ar
0
+ d
1
+ d
2
其中：上面的等式用到了万有引力定律和牛顿第二定律以及椭圆的相关知识。下面我们来确定近
月点和远月点的位置。
忽略月球自转。根据该阶段所处的状态，可以大体估计出其需要的时间约为 410s 左右。由
于嫦娥 3 号落点为 19.51
◦
W, 44.12
◦
N，根据假设，嫦娥三号软着陆轨道为过月球自转轴的平面，
且从南至北着陆，那么嫦娥三号的近月点的维度为 44.12
◦
− 12.6397
◦
= 31.4803
◦
，即嫦娥三号近
月点的经纬度为 (19.51
◦
W, 31.4803
◦
N)。根据对称性，远月点的经纬度为 (160.5
◦
E, 31.4803
◦
S)。
http://www.ma-xy.com 24 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
1.2.5 问题二的分析与求解
问题的分析
确定嫦娥三号的着陆轨道和在 6 个阶段的最优控制策略。
模型的建立与求解
(1) 主减速阶段
主减速阶段是距离月面 15m 到 3km，该阶段的主要任务是减速，在 3km 处速度降到 57m/s。
我们给出下面 4 种方案来求解此问题。
方案 1：2D 直角坐标系下的动力系统分析。以月心为坐标系原点，以月心指向近月点 (嫦
娥 3 制动的初始点) 方向为 y 轴方向，以垂直 oy 且指向娥 3 移动方向为 x 轴，建立平面直角坐
标系 xoy 如图 (1.15) 所示
图 1.15: 嫦娥 3 号软着陆平面直角坐标系
以 t 时刻为研究点进行受力分析，娥 3 主要受月球引力 G 和主发动机制动力 F 的影响，由
图 (1.15) 中的受力分析可知
G = gm
=
G
0
M
r
2
m
=
G
0
M
x
2
+ y
2
m
且有 α = arctan
x
y
。其中：G
0
= 6.672 × 10
−1
N 为万有引力常量，M = 7.350 × 10
22
kg 为月球质
量，m 为娥 3 的质量，且由于着陆过程中娥 3 消耗燃料制动，所以各时刻的娥 3 质量不同，记
为 m
t
，x, y 为坐标点、r 为月心距皆为时变量。
将娥 3 整体视为系统，我们主要关注系统在 t 时刻的系统变量：x
t
, y
t
, m
t
, F
t
, θ
t
, v
t
。为了便
于理解，不妨将时间离散化，研究系统下一时刻 (t + 1 时刻) 的系统值。如果我们已经知道了 t
http://www.ma-xy.com 25 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
时刻的系统值，那么 t + 1 时刻的系统值为
x
t+1
= x
t
+ (v
x,t+1
− v
x,t
)(t + 1 − t)
y
t+1
= y
t
+ (v
y,t+1
− v
y,t
)(t + 1 − t)
v
x,t+1
= v
x,t
+ a
x,t
(t + 1 − t)
v
y,t+1
= v
y,t
− a
y,t
(t + 1 − t)
a
x,t
=
F
t
sin θ
t
+ G
t
sin α
t
m
t
a
y,t
=
−F
t
cos θ
t
+ G
t
cos α
t
m
t
m
t+1
= m
t
−
F
t
v
e
α
t
= arctan
x
t
y
t
例如：我们给出初始值 x
0
, y
0
, v
0
, m
0
及 F
0
, θ
0
(待求变量)，即可得到系统下一时刻的状态值
x
t
, y
t
, v
t
, m
t
。不妨将上面的离散状态更新公式写为连续时间 t，有
˙x =
dx
dt
= v
x,t
˙y =
dy
dt
= v
y,t
˙v
x
=
dv
x
dt
= a
x,t
˙v
y
=
dv
y
dt
= a
y,t
a
x,t
=
F
t
sin θ
t
+ G
t
sin α
t
m
t
a
y,t
=
−F
t
cos θ
t
+ G
t
cos α
t
m
t
˙m =
dm
dt
= −
F
t
v
e
α
t
= arctan
x
t
y
t
http://www.ma-xy.com 26 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
由于速度是距离的导数，加速度是速度的导数，故上述系统方程可写为
¨x =
d
2
x
d
t
2
= a
x,t
¨y =
d
2
y
dt
2
= a
y,t
a
x,t
=
F
t
sin θ
t
+ G
t
sin α
t
m
t
a
y,t
=
−F
t
cos θ
t
+ G
t
cos α
t
m
t
α
t
= arctan
x
t
y
t
˙m =
dm
dt
= −
F
t
v
e
最优化控制问题：从连续时间 t 来看，我们需要给出推理 F
t
和夹角 θ
t
的函数，使系统满足
初始条件 (x
0
, y
0
, v
0
) 和终端条件 (x
t
, y
t
, v
t
)，并且使 t 时间消耗的能量最少；从离散时间 t 来看，
我们需要给出每时刻 t 的 F
t
, θ
t
，使系统满足初始条件 (x
0
, y
0
, v
0
) 和终端条件 (x
t
, y
t
, v
t
)，并且使
t 时间消耗的能量最少。上述内容可以写为如下数学模型 (为和一般文献相吻合，定义 0 时刻为
t
0
，终止时刻为 t
f
)
min
F,θ
J =
t
f
t
0
˙mdt = m
0
− m(t
f
)
s.t.
x(t
0
) = x
0
, y(t
0
) = y
0
, v
x
(t
0
) = v
x,0
, v
y
(t
0
) = v
y
0
, m(t
0
) = m
0
¨x =
d
2
x
dt
2
= a
x,t
¨y =
d
2
y
dt
2
=
a
y,t
a
x,t
=
F
t
sin θ
t
+ G
t
sin α
t
m
t
a
y,t
=
−F
t
cos θ
t
+ G
t
cos α
t
m
t
α
t
= arctan
x
t
y
t
˙m =
dm
dt
= −
F
t
v
e
x(t
f
) = x
f
, y(t
f
) = y
f
, v
x
(t
f
) = v
x,f
, v
y
(t
f
) = v
y,f
上述优化控制的可行性：如果是在连续时间 t 上讨论上述问题，则优化控制室一个 F
t
, θ
t
的
泛函问题，解 F
t
, θ
t
是无穷维的，很难求解，所以我们一般将时间 t 离散化处理。如果给定初始
状态 x
0
, y
0
, v
0
, m
0
，给定一个解 {F
t
}
t
f
t=t
0
, {θ
t
}
t
f
t=t
0
(不妨将解何为 u = [F
t
, θ
t
]
T
, {u
t
}
t
f
t=t
0
)，根据状
态更新公式 (即约束条件 s.t.)，我们可以求解 t 时刻的系统状态 {x
t
, y
t
, v
t
, m
t
}。并且要求在终端
t
f
时刻，系统状态满足要求 (x
f
, y
f
, v
f
)。如果有 {u
t
} 使系统达到要求 (s.t.)，我们即可得到最终
的目标值 J = m
0
− m
f
。
定义满足 s.t. 的解 {u
t
} 为可行解，其取值域为可行域 Φ，故而我们可以在可行域 Φ 中找到
{u
t
}
∗
，使得 J 最小。上面遗留的问题是：如何求解 {y
t
}
∗
。
http://www.ma-xy.com 27 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
方案 2：2D 极坐标系下的动力系统分析。在上面的直角坐标系的处理中，我们需要使 x, y
满足 t
f
时刻的月心距要求
r
f
=
x
2
f
+ y
2
f
并且，要将速度 v 在各时刻 t 上分解为 v
x
, v
y
。这样做虽然理解起来简单，但无疑饶了圈子，我
们不妨在极坐标下对系统进行动力学分析，因为在极坐标下，r 即为 r，v 与 r 方向垂直。构建
2D 极坐标系如图 (1.16) 所示
图 1.16: 嫦娥 3 号软着陆平面极坐标系
我们仍然假设¬不考虑其它天体对娥 3 的摄动；下降轨道与经纬线重合。我们所关心的系
统变量为 (r, θ, v, w, m, F, φ)，其中：F, φ 为系统控制变量；(r, v, θ, w, m) 为系统状态变量。
由运动分析，我们有
˙r =
dr
dt
= v
˙
θ =
dθ
d
t
= w
在切向上进行受力分析，飞行器受 F cos φ 及科里奥利力 的影响
F cos φ + 2vw = −m
dv
y
dt
其中：v
y
= wr。故有
˙w =
dw
dt
= −
1
r
(
F
m
cos φ + 2vw)
在径向 r 上进行受力分析
GMm
r
2
=
mv
2
r
F sin φ = m
dv
dt
其中：v = wr。故有
˙v =
dv
dt
=
F
m
sin φ −
GM
r
2
+ rw
2
http://www.ma-xy.com 28 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
由上述分析，我们得到娥 3 质心的运动方程为
˙r =
dr
dt
= v
˙
v
=
F
m
sin
φ
−
GM
r
2
+
rw
2
˙
θ =
dθ
dt
= w
˙w = −
1
r
(
F
m
cos φ + 2vw)
˙m = −
F
v
e
其中：v ∈ R 是娥 3 在矢径 r 方向上的速度；w ∈ R 是娥 3 在角 θ 上的角速度；m ∈ R
+
是娥 3
质量；G
0
为万有引力常量；M 为月球质量；v
e
为比冲。
设解为 u
t
= [F
t
, φ
t
]
T
，离散化后可写为 {u
t
}
t
f
t=t
0
，最优控制是寻找最优控制率 u
∗
t
，使着陆过
程消耗的燃料最小，并且满足初始条件和终端条件，即
min
u=F,φ
J =
t
f
t
0
˙mdt = m
0
− m(t
f
) = m
0
− m
f
s.t.
r
(
t
0
) =
r
0
, v
(
t
0
) =
v
0
, θ
(
t
0
) =
θ
0
, w
(
t
0
) =
w
0
, m
(
t
0
) =
m
0
˙r =
dr
dt
= v
˙
v
=
F
m
sin φ −
GM
r
2
+ rw
2
˙
θ =
dθ
dt
= w
˙w = −
1
r
(
F
m
cos φ + 2vw)
˙m = −
F
v
e
r(t
f
) = r
f
, v(t
f
) = v
f
, w(t
f
) = w
f
至此，主减速阶段的数学模型已经建立完毕。下面我们来讨论如何求解 {u
t
}
∗
。在此之前，先引
入科里奥利力，然后再介绍求解 {u
t
}
∗
的直接法和间接法。
科里奥利力：哥氏力又称科里奥利力，简称科氏力，是哥氏加速度的来源。哥氏加速度是由
于质点不仅作圆周运动，还做径向运动或轴向运动所产生的。当质量为 m 的质点相对于转动参
考系 (月球，角速度矢量 w) 的速度为 v 时，在转动参考系内观察到科里奥利力为
F
c
= 2mvw
科里奥利力只有在转动参考系中运动时才发生，科氏力方向垂直相对速度，不会改变相对速度大
小。
科里奥利力的简单推导：如图 (1.17) 所示
http://www.ma-xy.com 29 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
图 1.17: 科里奥利力示意图
质点 m 在转动参考系 S
′
中沿一光滑凹槽运动，速度为 v
′
，在惯性地面 S 中
F = m
(v
2
+ rw)
2
r
= m
v
2
r
+ 2mv
′
w + mrw
2
在非惯性系圆盘 S
′
中，向心加速度为
a
′
=
v
′2
r
w 的方向垂直 v 向上。
最优控制 u
∗
t
通常有 2 种方法：1. 直接法；2. 间接法。直接法是将 u
∗
t
离散化、参数化后进
行求解；间接法是在最优控制理论 (变分/极大值原理/动态规划) 的基础上，将优控问题变化为
终端 t
f
自由的两点边值问题，然后求解。
最优控制 u
∗
t
的直接法：由于月球软着陆的轨道优化的搜索空间是一个 u
t
∈ {u
t
} = R 函数
空间 (泛函空间)，而我们通常处理的最优化问题的搜索空间是一个 θ ∈ Θ 空间 (参数空间)。所
以我们无法用最优化数值方法直接得到轨道，这使得我们不得不另寻它法。既然我们需要参数函
数，不妨将 {u
t
} 设为参数函数，设
{
u
t
|
u
t
=
u
(
t
|
θ
)
, θ
∈
R
+
}
例如：我们可以设 u
t
为参数多项式、傅里叶多项式等各种含参函数形式
u
t
=
1
a + e
t
u
t
=
n
t=0
a
i
t
i
这样，我们不再是在函数空间 R 中寻找一个函数 u
∗
t
使 J 最优，而是在参数空间 Θ 中寻找一个
参数 θ，进而有 u
∗
t
使 J 最优。一般文献中都使用含参多项式的形式作为 u
t
的替代，我们下面也
采用这种形式。设 u
t
是一高 n 次多项式
u
t
=
n
i=0
a
i
t
i
其中：{a
i
}
n
i=0
= θ ∈ Θ，共引入 n + 1 个参数。
http://www.ma-xy.com 30 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
在引入 u
t
之后，再要考虑的是：将 u
t
以连续函数的形式带入到动力方程当中进行求解，还
是将时间 t 离散化逐点求解。我们先考虑将时间 t 离散化处理，将 [t
0
, t
f
] 离散化为 N 段，由于
t
f
不定，设
t
i
= t
0
+
i
N
(t
f
− t
0
)
其中：N 为时间段，时间段 [t
i
, t
i+1
], i = 1, 2, . . . , N 的两端称为节点。这样，我们就有 i 段 [t
i−1
, t
i
]
对应的 [u
t
i−1
, u
t
i
]
(
{
a
i
}
给定后，
u
t
对应具体数值
)
，进而可以得到
[
x
t
i−1
, x
t
i
]，进而得到 x
t
f
以及
J。我们求最优 {a}, t
f
使 J 最优。
上述方法可以用智能算法 GA 等进行求解，当然，其它优化算法 BFGS 等亦可。而且注意到，
当 F
t
为固定值时，目标变为 min J = t
f
。在具体求解过程中，我们需要给出参数的初始值/估计
值/迭代出发点。t
f
的估计还简单一些 (可以自行建模估计)，但到了 {a
i
} 却没有实际的物理意
义，这使得他们的估计量是“天马行空”的。如果系统 x
t
对参数 {a
i
} 是及其敏感的，那么，{a
i
}
初始估计将成为问题所在。
问：如何解决 {a
i
} 初始估计无依据 (无物理意义)？缺：未讨论 t 不离散的情况。
插值-分段三次埃米尔特 Hermite-Simpsorn。上面，我们将时间 t 离散后，设 u
t
为多项式，
那么，我们能不能假设 x
t
是多项式。承继上面的时间离散化，设 h
i
= t
i
− t
i−1
, i = 1, 2, . . . , N，
我们并不打算在 [t
0
, t
f
] 内设 x
t
=
ˆ
f(t)，而是在某一小时间段 [t
i−1
, t
i
] 进行设置，设
x = a
0
+ a
1
t + a
2
t
2
a
3
t
3
=
ˆ
f(t|θ)
其中：t ∈ [t
i−1
, t
i
]，a 承继上面为参数 θ。
也就是说，x 在 [t
i−1
, t
i
] 内皆是三次多项式，但不同时段的参数 θ 的取值不同，即 x 为分段
三次多项式。不失为一般性，我们在时段 [t
0
, t
1
] 内研究 x。
x = a
0
+ a
1
t + a
2
t
2
a
3
t
3
(1.3)
t ∈ [t
0
, t
1
]
假设我们知道了以下条件
x(t
0
) = x
0
x(t
1
) = x
1
x
′
(t
0
) = x
′
0
= f(x
t
0
)
x
′
(t
1
) = x
′
1
= f(x
t
1
) (1.4)
其中：˙x = f (x
t
, u
t
, t)。当 x
0
给定后，只要有 u
0
即有 x
′
0
。那么，我们就可以在 [t
0
, t
1
] 区间内对
函数 x(t) 进行插值，例如用三次 Hermite 多项式插值 ([t
0
, t
f
] 段上即为分段三次 Hermite 多项
式插值)。
注：关于插值的数值分析技术，可以参考《数值分析》李乃成 P127.
http://www.ma-xy.com 31 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
式 (1.3) 有四个参数待求，而我们有 4 个已知条件 (1.4)，我们可以利用解方程组的形式求解
参数 a
i
x
1
˙x
1
x
2
˙x
2
=
1 0 0 0
0 1 0 0
1 1 1 1
0 1 2 3
a
0
a
1
a
2
a
3
解上述方程组即可得到 {a
i
}
3
i=0
，将 a
i
带入到 H 三次多项式有
H
3
(t) = ˆx(t) =
1 + 2
t
0
− t
t
0
− t
1
t
−
t
1
t
0
− t
1
2
+
1 + 2
t
1
−
t
t
1
− t
0
t
−
t
0
t
1
− t
0
2
x(t
1
)
+ (t − t
0
)
t − t
1
t
0
− t
1
2
x
′
(t
0
) + (t − t
1
)
t − t
0
t
1
− t
0
2
x
′
(t
1
)
其中：t
1
− t
0
= h
1
， ˙x = f (x, u, t)。
这样，我们就得到了 [t
0
, t
1
] 区间内 x(t) 的 Hermite 插值估计 H
3
(t)。但是，由于我们需要
求解 x
0
, x
1
, u
0
, u
1
，每给一组 (x
0
, x
1
, u
0
, u
1
)，都会有一个 H
3
(t) 与之对应，x
0
, x
1
, u
0
, u
1
取不同
值时的 H
3
(t) 示意图如图 (1.18) 所示
图 1.18: x-0 等取不同值时的 H-3 示意图
既然有这么多 H
3
(t)，那么我们需要约束 H
3
(t)，不能让其乱跑。我们用条件 ˙x = f(x, u, t)
来进行约束。任给 t ∈ [t
0
, t
1
]，u(t) 给定后，如果知道 x(t)，则 ˙x(t) 是知道的，H
3
(t), H
′
3
(t) 是
已知的。我们可以在 [t
0
, t
1
] 内取 k 个点 {t
k
}，要求这 k 个点的 H
′
3
(t) 和 ˙x(t) = f(H
3
(t), u(t), t)
相等或者近似 (配点)。为计算简单，这里 k = 1，置 t
k
=
t
0
+t
1
2
(当 k 越大点数越多时，计算越精
确，但计算量会相应变大)
H
3
(t
k
) =
x
0
+ x
1
2
+
h
1
8
(f
0
− f
1
)
H
′
3
(t
k
) = −
3(x
0
− x
1
)
2h
1
−
f
0
− f
1
4
其中： ˙x
0
= f(x
0
, u
0
, t
0
) = f
0
， ˙x
1
= f(x
1
, u
1
, t
1
) = f
1
。
在给出 u
k
后，要求 H
′
3
(t
k
) = ˙x(t
k
)，这里我们两种方案。方案 (1)，我们利用 H
3
(t
k
) 直接
替代 x(t
k
)，有
H
′
3
(t
k
) = ˙x(t
k
) = f(H
3
(t
k
), u(t
k
), t
k
)
方案 (2)，三阶 Simpson 方法然后再用 H
3
(t
k
) 直接替代 x(t
k
)，区间 [a, b] 上的三阶 Simpson 求
积公式如下
S
i
=
b − a
6
[f(a) + 4f (
a + b
4
) + f(b)]
http://www.ma-xy.com 32 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
结合上面的三阶 Simpson 求积公式，我们有
x
i+1
= x
i
+
h
i
6
[f(x
i
, u
i
, t) + 4f (H
3
(t
k
), u(t
k
), t
k
) + f(x
i+1
, u
i+1
, t
i+1
)]
x
1
= x
0
+
h
1
6
[f(x
0
, u
0
, t
0
) + 4f(H
3
(t
k
), u(t
k
), t
k
) + f(x
1
, u
1
, t
1
)]
注：x 为 f(t)，H
0
= x
0
= x(0), H
1
= x
1
= x(1) 一定是高度重合的。
经过上面的分析，在 [t
0
, t
1
] 时段内的处理我们已经基本完成，很容易将上述结论推广到 N
段。在 [t
0
, t
1
] 段内，我们假设已知 (x
0
, x
1
, ˙x
0
, ˙x
1
, u
t
k
)，故这样变为我们需要求解的优化变量，由
于 ˙x 与 u 有关，故我们需要求解的优化变量为
z = [x
T
0
, u
T
c0
, u
T
0
, x
T
1
, u
T
c1
, u
T
1
, . . . , x
T
N
, u
T
cN
, u
T
N
]
其中：u
T
ci
即为第 [t
i−1
, t
i
] 段内的配点 t
k
的控制量。
这种方法将最优控制问题离散化、参数化后，转化为非线性规划问题，如下
min
z
J =
t
f
t
0
˙mdt = −m(t
f
)
s.t.
r(0) = r
0
, v(0) = v
0
, θ(0) = 0, w(0) = w
0
, m(0) = m
0
u
min
⩽ u ⩽ u
max
˙x = f (x, u, t)
r(t
f
) = r
f
, v(t
f
) = v
f
, w(t
f
) = 0
上述 s.t. 并不易于书写，我们来捋一下：¬首先，哪些是待优化变量？
z = [x
T
0
, u
T
c0
, u
T
0
, x
T
1
, u
T
c1
, u
T
1
, . . . , x
T
N
, u
T
cN
, u
T
N
]
将 [t
0
, t
f
] 等分为 N = 20 段，共 21 个节点。各节点上有 5 个状态变量 x，2 个控制变量 u。20
段中的中间配点有 2 个 u(k = 1)，故待优化变量个数为
(N + 1)(n
x
+ n
u
) + N(n
u
)
其中：n
x
为 x 变量个数；n
u
为 u 变量个数，N 为段数。等式约束。上面的优化模型中的等式
约束是
初始条件 + 终端条件 + (Nn
x
)的 HS 等式约束
®不等式约束为
u
min
⩽ (N + 1)n
u
⩽ u
max
u
min
⩽ (N)n
u
⩽ u
max
承继上面的想法，能否直接给出 (x
i
, x
i+1
, u
i
, u
i+1
)，并设置一定的等式条件。文
4利用 4 阶
R-K(Runge-Kutta) 法实现了直接法。这里我们给出 R-K 公式 (可以参考3P273、2P9、1P141)，
http://www.ma-xy.com 33 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
考虑如下一阶常微分初值问题
y
′
(
x
) =
f
(
x, y
(
x
))
y(a) = y
0
上述问题的 4 阶 R-K 计算公式为
y
i+1
= y
i
+
1
6
(K
1
+ 2K
2
+ 2K
3
+ K
4
)
K
1
= hf(x
i
, y
i
)
K
2
= hf(x
i
+
1
2
h, y
i
+
1
2
K
1
)
K
3
= hf(x
i
+
1
2
h, y
i
+
1
2
K
2
)
K
4
= hf(x
i
+ h, y
i
+ K
3
)
其中：h 为 x 的步长。改进的 4 阶 R-K 公式为
y
i+1
= y
i
+
1
6
(K
1
+ 4K
3
+ K
4
)
K
1
= hf(x
i
, y
i
)
K
2
= hf(x
i
+
1
2
h, y
i
+
1
2
K
1
)
K
3
= hf(x
i
+
1
2
h, y
i
+
1
4
K
1
+
1
4
K
2
)
K
4
= hf(x
i
+ h, y
i
− K
2
− 2K
3
)
以及
y
i+1
= y
i
+
1
8
(K
1
+ 3K
2
+ 3K
3
+ K
4
)
K
1
= hf(x
i
, y
i
)
K
2
= hf(x
i
+
1
3
h, y
i
+
1
3
K
1
)
K
3
= hf(x
i
+
2
3
h, y
i
−
1
3
K
1
+ K
2
)
K
4
= hf(x
i
+ h, y
i
+ K
1
− K
2
+ K
3
)
承继上面的 [t
0
, t
f
] 离散化方法，并且有系统状态方程
˙x = f (x, u, t)
在区间 [t
i
1
, t
i
], i = 1, 2, . . . , N (t
N
= t
f
)，假设我们有 (x
i−1
, u
i−1
, x
i
, u
i
)，那么由 4 阶 R-K，我们
有 x
i−1
, x
i
的计算公式
x
i
=
x
i−1
+
1
6
(
K
1
+ 2
K
2
+ 2
K
2
+
K
4
)
http://www.ma-xy.com 34 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
其中：
K
1
= h
k
f(x
k
, u
k
, t
k
)
K
2
= h
k
f(x
i
+
1
2
K
1
,
¯
λ
k+1
,
¯
t)
K
3
= h
k
f(x
i
+
1
2
K
2
,
¯
λ
k+1
,
¯
t)
K
4
= h
k
f(x
i
+ K
3
,
¯
λ
k+1
, t
k+1
)
h
k
= t
k+1
− t
k
, h
i−1
= t
i
− t
i−1
, i = 1, 2, . . . , N
¯
t =
1
2
(t
k
+ t
k+1
)
λ
k+1
= λ(
¯
t)
经过上述离散化处理后，原优化控制问题变为非线性规划问题，其优化变量为
z = [x
0
, u
0
, x
1
, u
1
, . . . , x
N
, u
N
]
T
上面两种方法都将最优控制问题处理为含约束的非线性规划问题，处理非线性规划问题的方
法有很多，可以参考前面相应的章节。非线性规划的处理方法一般有：1. 可行方向法；2. 罚函数
法 (内点法)；3. 梯度投影法；4. 逐步二次规划 SQP。而 MATLAB 主要采用 SQP 方法。
最优控制 u
∗
t
的间接法：跋山涉水，终于到了最优控制 u
∗
t
的间接解法 (后面，我们还会建立
软着陆的
3D
动力系统
)
。前面，我们花了很大篇幅来讨论系统
˙
x
=
f
(
x, u, t
)
，泛函
J
(
x
t
, u
t
, t
)
。
泛函极值存在条件是欧拉方程、最优控制中的控制量 u
t
无约束时的变分法、控制量 u
t
有约束时
的 P 极大值原理及动态规划。下面，我们就要用这些理论来分析娥 3 轨道的优化问题。
我们可以看出，娥 3 轨道最优控制问题是一个控制量 u = (F, φ) 有约束且约束 t
f
自由、
x(t
f
) 有约束的最优控制问题，因此，我们可以用 P 极大值原理或动态规划进行求解。下面，我
们将用 P 极大值原理求解最优控制 u
∗
= (F, φ)
∗
，利用 P 极大值原理求解最优控制律 u
∗
t
，从而
将问题转化为两点边值问题，当然也可以进行离散化求解。
1962.Kriegsman 和 Reiss 利用变分法求解最优控制律 u
t
。2000 年王大秩等利用 Pontryagin
极大值原理求解 u
t
，将问题转变为两点边值问题，并用一种初值猜测的打靶算法进行求解。2004
年，单永正等利用 Pontryagin 极大值原理，将问题转化为两点边值问题，并将两点边值问题转
化为非线性规划问题进行求解。2004 年，赵松吉等利用初值猜测技术和线性摄动求解两点边值
问题。
下面，我们用 Pontryagin 极大值原理求解最优控制律 u
t
= (F
t
, φ
t
)
T
。记
f(x, u, t) =
f
1
= v
f
2
=
F
m
sin φ −
G
0
M
r
2
+ rw
2
f
3
= w
f
4
= −
1
r
(
F
m
cos φ + 2vw)
f
5
= −
F
v
e
http://www.ma-xy.com 35 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
于是动力系统方程可以写为
˙x = f (x, u, t) = f (x(t), u(t), t)
(1) 引入拉格朗日乘子/共轭变量 λ = [λ
r
, λ
v
, λ
θ
, λ
w
, λ
m
]
T
，构建 H(哈密顿函数)
H(x, λ, u, t) = λ
T
f(x, u, t)
= λ
r
f
1
+ λ
v
f
2
+ λ
θ
f
3
+ λ
w
f
4
+ λ
m
f
5
=
F
mr
rλ
v
sin φ − λ
w
cos φ −
λ
m
v
e
mr
+ H
2
= H
1
(F, φ) + H
2
其中：H
2
表示与控制变量 F, φ 无关部分。由 P 极大值原理，有 max H ≜ max H
1
。
(2) 写正则方程，有
˙
λ = −
∂H
∂x
(= 0)
˙x = −
∂H
∂λ
(= 0) = f(x, u, t)
推得
˙
λ
r
= −
2λ
v
G
0
M
r
3
− λ
v
w
2
−
λ
w
F cos φ
mr
2
− 2λ
w
vw
˙
λ
v
=
2wλ
w
r
− λ
r
˙
λ
θ
= 0
˙
λ
w
=
2vλ
w
r
− λ
θ
− 2rwλ
v
˙
λ
m
=
λ
v
F sin φ
m
2
−
λ
w
F cos φ
rm
2
(3) 写控制方程，有
˙u =
∂H
∂u
(= 0)
(4) 终端约束
G =
r(t
f
) = r
f
v(t
f
) = v
f
w(t
f
) = w
f
= 0
(5) 横截条件
λ(t
f
) =
∂J
∂x
+
∂G
∂x
r
t=t
f
其中：r 为拉格朗日乘子。
http://www.ma-xy.com 36 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
(6) 控制律约束
F
min
⩽ F
i
⩽ F
max
φ
min
⩽ φ
i
⩽ φ
max
由于求 u 上的极值，令控制方程等于 0
˙u =
∂H
∂u
= 0
有 φ, F 的最优控制律
φ
∗
= − arctan
λ
v
r
λ
w
F
∗
=
F
max
s
(
t
)
>
0
0 s(t) < 0
其中：s(t) =
λ
v
sin φ
m
−
λ
w
cos φ
mr
−
λ
m
v
e
。
在求解最优控制律 u
∗
t
后，要求 x
∗
, λ
∗
使 H 最大。将 u
∗
t
带入正则方程当中，利用初始条件
和终端条件对正则方程积分，即可求出最优 x
∗
, λ
∗
，此时，问题转化为微分方程中的两点边值问
题。
简单分析一下，在离散化下考虑，在 x(t
0
) 已知的情况下，给定一个初始值 λ(t
0
)，我们会得
到相应的最优控制律 u(t
0
)
∗
，进而由 ˙x = f 得到 x(t
1
)。依次迭代，可以得到终端 x(t
f
) = x
f
(在
给定初始 t
f
后)。进一步观察 λ
m
的取值，对 λ
/m
及 x 没有影响，因此，只要先选取迭代初始值
[λ
r
, λ
v
, λ
w
]
t=t0
，确定 u
t
0
，再给出已知 x(t
0
) = x
0
。在正则方程上积分，即可得到末端状态 x(t
f
)，
即 x(t
f
) 是初值 [λ
r
, λ
v
, λ
w
]
t=t0
, t
f
的函数，进而性能 J 也是是 (λ
r
(t0), λ
v
(t0), λ
w
(t0), t
f
) 的函
数。这样，两点边值问题就被离散化为非线性规划 (只有等式约束)。
迭代初始值 t
f
可以根据一定的公式进行估计，文5给出估计
t
f
=
∆m
| ˙m|
=
∆m
F /v
e
=
m
0
[1 − exp(−∆v/v
e
g)]
F /v
e
=
m
0
[1 − exp(−
r
2
0
+ 2gh/v
e
g)]
F /v
e
其中：| ˙m| 为发动机燃料秒耗量。
但是迭代初始值 (λ
r
(t0), λ
v
(t0), λ
w
(t0), t
f
) 的选取则没有那么多依据，我们只能凭空设置。
如果正则方程对初始值仍然敏感，那么 (λ
r
(t0), λ
v
(t0), λ
w
(t0)) 将成为问题的所在。下面的“伴
随-控制”变换可以解决这一问题。设
λ
v
(t
0
)r
0
= sin(φ
0
)
λ
w
(t
0
) = − cos(φ
0
)
http://www.ma-xy.com 37 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
将 tan φ = −
λ
v
r
λ
w
对
t
求导，有
˙φ
cos
2
φ
= −
˙
λ
v
r
λ
w
+
λ
v
r
λ
w
−
λ
v
r
˙
λ
w
λ
2
w
(1.5)
将正则方程和 λ
0
v
, λ
w
(t
0
) 带入式 (1.5)，有
λ
r
=
˙φ
λ
w
r
+
2λ
w
w
r
−
λ
v
v
r
+
2λ
v
rw
λ
w
于是有
λ
r
(t
0
) =
˙
λ
λ
w
(t
0
)r
0
+
2λ
w
(t
0
)w
0
r
0
−
λ
v
(t
0
)v
0
r
0
+
2λ
v
(t
0
)r
0
w
0
λ
w
(t
0
)
这样，(λ
r
(t0), λ
v
(t0), λ
w
(t0)) 就可以用 φ
0
, ˙φ
0
替代，可以取合理的初始值了。
引：两点边值问题的数值解。考虑下面一个简单的二阶常微分方程的两点边值问题
y
′′
(x) = f(x, y, y
′
)
y(a) = y
a
y(b) = y
b
a < x < b
前面我们讨论了常微的初值问题，我们能不能在初值问题的基础上求解边值问题？能，但对于二
阶常微的初值 y(a) = ∆
1
, y
′
(a) = ∆
2
，我们还不足以求解边值问题，我们缺少一个条件。那我们
就设一个条件好了，我们设 y
′
(z) − m，并舍去 y(b) = y)b，有
y
′′
(x) = f(x, y, y
′
)
y(a) = y
a
y
′′
(a) = m
数值求解上述两点边值问题。一个正常的思路是，我们不断调整 m，使 y(b) = y
b
，如图 (1.19)
所示
图 1.19: 两点边值问题的打靶方法
两点边值问题的打靶方法算法伪代码如 (1) 所示
http://www.ma-xy.com 38 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
算法 1 2 点边值打靶法 ode
1: 初始化：m，ε，t
max
，i := 0。
2: while t = 1, 2, . . . , t
max
do
3: // 解初值问题
y
′′
(x) = f(x, y, y
′
)
y(a) = y
a
y
′′
(a) = m
i
4: β
i
= y(b)。
5: if |β
i
− β| > ε then
6: 修正 m
i+1
。
7: end if
8: end while
上述伪代码中遗留的问题是：如何修正 m。理想的 m 是 y(b, m) = y
b
，但在数值求解过程
中，不存在解析表达式，并不易解。一种可行的方式是取
m
2
=
y
b
y
b1
m
1
此后，由非线性方程求根的割线法更新 m
i
m
i
= m
i−1
−
y(b, m
i−1
) − y
b
y(b, m
i−1
) − y(b, m
i−2
)
(m
i−1
− m
i−2
) i = 3, 4, . . .
Matlab 中的 ode45 用于求解初值问题，bvp4c 用于求解边值问题。注：由于状态变量的量
级差异大，在轨道积分过程中会导致有效位数的损失，归一化处理可以克服这一缺点。在一般论
文中皆有归一化说明，这里就不做叙述了。
方案 3：3D 球面坐标系下的动力系统分析。很明显的一个问题是：我们用上面的方案 1 或
者方案 2 中的 2D 动力学模型只能做某一方向上的平移，而不能做平面上的移动。而题中既然给
了空间上的数字高程图，那我们就有必要在 3 维空间中建立着陆轨道。当然，我们可以在未拍照
之前的主减速阶段和快速调整阶段建立 2D 模型，在拍照之后建立 3D 模型，但这无疑增加了复
杂性。幸运的是，3D 模型的建立和求解与 2D 模型是相似的。
如前面 2D 分析的那样，我们要考虑的问题是：1. 在直角坐标系下建模还是在球面坐标系下
建模；2. 是否考虑月球自转；3. 用直接法还是间接法求解最优控制问题。3D 球面坐标系下的动
力学模型如图 (1.20) 所示
http://www.ma-xy.com 39 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
图 1.20: 3D 球面坐标系下的动力学模型
假设：1. 忽略月球自转；2. 月球引力场为均匀引力场，g 为常值；3. 不考虑其他星体对娥 3 的
摄动。以月心 o
1
为坐标原点，
−−→
o
1
x
1
指向近月点 (娥 3 下降点)，
−−→
o
1
y
1
指向娥 3 着陆方向，
−−→
o
1
z
1
按
右手法则确定，建立月心惯性坐标系 o
1
x
1
y
1
z
1
。o
0
o
1
为娥 3 到月心的距离 r，α 为 o
0
o
1
在 y
1
o
1
z
1
面投影与
−−→
o
1
y
1
正方向夹角，β 为 o
0
o
1
与
−−→
o
1
x
1
正方向夹角，由此有球面坐标系 (r, α, β)。
此外，我们还需要确定推力 F 的方向。以娥 3 质心 o
0
为原点，
−−→
o
0
x
0
为背月心
−−→
o
1
o
0
延长线，
−−→
o
0
y
0
指向娥 3 运行方向，
−−→
o
0
z
0
由右手法则确定。建立着陆轨道坐标系 o
0
x
0
y
0
z
0
，ψ, φ 为推力 F
的方向角，形成球面坐标系 (F, ψ, φ)。
在上面的假设下，月球软着陆的动力学方程可以表示为
¨
a =
¨
r =
F
m
−
r
r
3
G
0
M
其中：
¨
r 是娥 3 的加速度矢量，
F 为发动机动力矢量，r 为娥 3 的位移矢量，G
0
为万有引力常
量，M 为月球质量。不妨令 G
0
M = µ(常量)。
将位移矢量 r 按 o
0
x
0
y
0
z
0
方向轴分解为 r
x
, r
y
, r
z
；将 r 按 o
1
x
1
y
1
z
1
方向轴分解为 r
′
x
, r
′
y
, r
′
z
。
通过 o
0
x
0
y
0
z
0
我们可以看到，o
0
x
0
y
0
z
0
沿 o
0
x
0
轴逆方向转动 (90 − α)
◦
，然后沿 o
0
y
0
轴逆时针
转动 β，最后沿 o
0
z
0
轴顺时针转动 (90 − α)
◦
即可与 0
1
x
1
y
1
z
1
归正。因此，从着陆轨道坐标系
o
0
x
0
y
0
z
0
转到月心惯性坐标系 o
1
x
1
y
1
z
1
的转动矩阵 I 为
I =
cos α cos β sin β − sin α cos β
− cos α sin β cos β sin α sin β
sin α 0 cos α
即 [r
′
x
, r
′
y
, r
′
z
]
T
= I[r
x
, r
y
, r
z
]
T
。
http://www.ma-xy.com 40 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
对娥 3 质心 o
0
沿 o
1
x
1
y
1
z
1
各方向轴进行受力分析，用 a
z
, a
y
, a
z
表示各方向轴的加速度，有
¨r − ˙α
2
r sin
2
β −
˙
β
2
r = a
x
2
˙
β ˙r +
¨
βr − ˙α
2
r cos β sin β = a
y
2 ˙α ˙r sin β + ¨αr sin β + 2 ˙α
˙
βr cos β = a
z
a
x
=
F cos ψ
m
−
µ
r
2
a
y
=
F sin ψ cos φ
m
a
z
=
F sin ψ sin φ
m
˙m = −
F
v
e
用 v
x
, v
y
, v
z
表示速度 v 在 o
1
x
1
y
1
z
1
各方向轴的分量，化简上述方程，有
˙r = v
x
˙
β =
v
y
r
˙α =
v
z
r sin β
˙v
x
=
F cos ψ
m
−
µ
r
2
+
v
2
y
+ v
2
z
r
˙v
y
=
F sin ψ cos φ
m
−
v
x
v
y
r
+
v
2
z
r tan β
˙v
z
=
F sin ψ sin φ
m
−
v
x
v
y
r
−
v
y
v
z
r tan β
˙m = −
F
v
e
其中：µ = G
0
为常量，v
e
为发动机比冲常量。由此，我们得到了 3D 球面坐标系下的轨道运动
方程 (系统)。
方案 4：3D 直角坐标系下的动力系统分析。在直角坐标系下对系统进行受力分析是简单的。
设 v
x
, v
y
, v
z
沿用上面的设置，则有
˙x = v
x
˙y = v
y
˙
z
=
v
z
设 v
′
x
, v
′
y
, v
′
z
为速度 v 在 o
0
x
0
y
0
z
0
各方向轴分量，则有
˙v
′
x
=
F cos ψ
m
−
µ
r
2
˙v
′
y
=
F sin ψ cos φ
m
˙v
′
z
=
F sin ψ sin φ
m
http://www.ma-xy.com 41 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
同时，由坐标转换有
[ ˙v
x
, ˙v
y
, ˙v
z
] = I[ ˙v
′
x
, ˙v
′
y
, ˙v
′
z
]
于是，3D 直角坐标系下的动力方程 (系统) 为
˙x = v
x
˙y = v
y
˙z = v
z
˙v
′
x
=
F cos ψ
m
−
µ
r
2
˙v
′
y
=
F sin ψ cos φ
m
˙v
′
z
=
F sin ψ sin φ
m
[ ˙v
x
, ˙v
y
, ˙v
z
] = I[ ˙v
′
x
, ˙v
′
y
, ˙v
′
z
]
˙m = −
F
v
e
文献6建立了 3D 直角坐标系且考虑了月球自转的着陆动力系统。在上面不考虑月球自转的
前提下，我们建立了软着陆的 3D 直角坐标系和球面坐标系下的动力系统。下面，我们给出此动
力系统的最优控制模型。
设 x = [r, β, α, v
x
, v
y
, v
z
, m]
T
为系统的状态变量，u = [F, ψ, φ] 为系统的控制变量，设 ˙x =
f(x, u, t) 为系统的状态方程，其中 f 为标量函数。我们的目标是控制 u，使系统消耗燃料最少，
故娥 3 的最优控制模型为
min
u
J =
t
f
t
0
˙mdt = m
0
− m
f
s.t.
r(t
0
) = r
0
, β(t
0
) = β
0
, α(t
0
) = α
0
,
v
x
(t
0
) = v
x0
, v
y
(t
0
) = v
y0
, v
z
(t
0
) = v
z0
, m(t
0
) = m
0
˙x = f (x, u, t)
r(t
f
) = r
f
, β(t
f
) = β
f
, α(t
f
) = α
f
, v
x
(t
f
) = v
xf
, v
y
(t
f
) = v
yf
, v
z
(t
f
) = v
zf
F
min
⩽ F < F
max
φ
min
⩽ φ ⩽ φ
max
ψ
min
⩽
ψ
⩽
ψ
max
其中：J 为积分型指标泛函。
上述最优控制问题是一个终端给定、t
f
不定、系统状态含约束、控制函数 u 含约束的积分
型最优控制问题。像分析 2D 最优控制那样，其解法可分为直接法和间接法。(1) 直接法：我们可
以将 [t
0
, t
f
] 离散化后，将泛函最优控制问题转化为参数最优化问题进行求解。(2) 间接法：利用
P 极大值原理，将最优控制问题转化为两点边值问题，继而可以用打靶算法和离散化参数化优化
http://www.ma-xy.com 42 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
算法进行求解。哈工大. 张仲满 2009.6 硕士论文中提到用 PID 自适应网络求解上述最优控制问
题，同时讨论了 P 极大值原理 + 殆尽的打靶算法在该问题中的应用；哈工大. 单永正.2009.7 博
士论文中提到用参数化技术 + 约束变换额方法来参数化求解最优控制问题，文中给出了相应的
理论与算法，并用 miser3.2 最优控制软件进行求解；哈工大. 张锋.2009.6 硕士论文中分析了姿态
发动机的姿态控制，并用 Lyapunov 分析系统稳定性，提出在软着陆末端进行姿态发动机和主减
速发动机联合控制策略；关于模糊、灰色、随机系统的最优控制问题可以参考葛宝明, 林飞《先
进控制理论》。由于 3D 动力系统的求解和 2D 系统相似，这里，我们不再继续 3D 求解工作。
(2) 快速调整阶段
(3) 粗避障阶段
在减速阶段和快速调整阶段之后，娥 3 进入粗避障和精避障阶段。粗避障阶段的范围是月面
2.4km 到 100m 区间，其主要目的是避开大的陨石坑，实现在设计着陆点上方 100m 处悬停，并
初步确定落月地点。嫦娥 3 号在月面 2.4km 时，对正下方 2300 × 2300m 的范围进行拍照，获得
数字高程地图如图 (1.7)。
下面，我们将利用数字高程图中的信息来建立避障模型。由于拍摄区域相对月球很小，不妨
假设地面为水平，不考虑月球弧度。以娥 3 所在位置俯视月面，进行拍摄，以其右下方顶点为 o
点，以高程图 0 高度水平面为 xoy 面，建立直角坐标系如图 (1.21) 所示。其中：(x
0
, y
0
) 为嫦娥
3 在月面正下方投影，图中阴影表示陨石坑，圆叉表示山峰。由于是数字高程图，我们可以得到
2300 × 2300 中各点的高度，不妨设各点高度为 h(x
i
, y
j
), i, j = 1, 2, . . . , n，n 为行像素点个数。
图 1.21: 避障坐标系
由于地面上有许多噪声，我们不妨在进行选择区域之间对图像进行去噪处理。之后，我们考
虑着陆点的条件，我们应该综合这些信息来评价某一点 (x
i
, y
j
)。我们考虑以下 3 个评价指标：
1. 点 (x
i
, y
j
) 应该平缓；
2. 点 (x
i
, y
j
) 平缓且其所在地区也具有一定的平缓度；
3. 点 (x
i
, y
j
) 应该距离 (x
0
, y
0
) 不远，以避免消耗较多燃料。
我们在 (x
i
, y
j
) 中找到最优着陆点为 (x
1
, y
1
)。我们有了上述 3 个评价之后，接下来的工作就
是量化指标，并构建评价函数。我们应该逐点 (x
i
, y
j
) 来考虑问题，即给出每个点的评价值，并
从中选优，而不是划大的区域，那样我们找不到一个确定的点来指导娥 3 着陆。
http://www.ma-xy.com 43 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
(一) 着陆区平缓性指标。首先，(x
i
, y
j
) 只有在靠近水平面上时，才平缓一点，也就是说，高
度 h
i
较大较小都不好。我们可以用 h
i
的分布来看一下 h
i
的情况。其次，与平缓最接近的数学
定义是梯度。我们用 Sobel 算子计算高程图的梯度函数 S(x, y)，并进行降噪处理，颜色越深地区
越平坦。定义平稳性因子为
P (x, y) = |h| + S(x, y)
P (x, y) 越小越好。
(二) 着陆区耗能型指标。点 (x
i
, y
j
) 距离 (x
0
, y
0
) 越远，则所耗能量越多。不妨设置耗能函
数为
F (x, y) = (x − x
0
)
2
+ (y − y
0
)
2
(三) 着陆区广阔性指标。由于飞行器娥 3 具有一定的大小/体积，并且在控制下落过程中，有许
多误差的干扰，会导致着陆点并不是很准确，所以我们需要着陆点具有一定的广阔性。以着陆点
(x
i
, y
j
) 为圆心，r
i
为半径的区域内的平均平缓都要在
D(x, y) = E(p|r)
对于一个地点 (x
i
, y
j
) 的评价要综合上述 3 点。我们设置 3 个指标的权重为 λ
1
, λ
2
, λ
3
，则最
终个的评价指标函数为
A(x, y) = λ
1
P (x, y) + λ
2
F (x, y) + λD(x, y|r)
(4) 精避障阶段
精避障阶段的区间是月面 100m 到 30m，要求娥 3 悬停在 100m 处时对着陆区 100m 范围
内进行拍摄，并获取三维数字高程图，如图 (1.8) 所示。分析三维数字高程图，避开大的陨石坑，
确定最终着陆点，实现在着陆上方 30m 处水平方向速度为 0m/s。
我们介绍一种基于李雅普诺夫稳定控制理论的非线性方法。李雅普诺夫控制制导法主要利用
李雅普诺夫稳定理论来进行控制律 u 的求解。下面介绍李雅普诺夫稳定判定定律
定理 (李雅普诺夫稳定判定定律) 设系统状态方程为
˙x = f (x, u, t)
其中：f(0, 0, t) = 0。如果存在一具有连续一阶偏导数的标量函数 (即李雅普诺夫函数)v(x, u, t)，
并且满足
v(x, u, t)为正定
˙v(x, u, t)为负定
则状态空间原点处的平衡状态是一致渐进稳定的。
http://www.ma-xy.com 44 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
选取正定的李氏函数 v，要求李氏函数 v 即能表示当前位置危险程度，又能表示距离着陆点
的位置关系。这样，v 函数就代表了飞行器娥 3 的期望性能。令他的导数为负，求取控制律，可
以保证飞行器的状态达到期望位置 (即使系统在期望状态处是一只渐进稳定的)。在月球精避障
阶段 (100m-30m)，控制律的目的是保证娥 3 在指定的高度 (30m) 处平移到安全着陆点上方。飞
行器平移的初始点为悬停点 (100m)，其末端点位于安全着陆点上方 (30m)，并且希望飞行器到
达着陆点上方时的平移速度为 0。
要求正定的李氏函数 v，即能表示当前状态/位置危险程度，又能表示到预定着陆点的燃料
消耗情况，故构建由能量函数和势危函数组成的 v，如下：¬能量函数
ϕ
p
(x, y, z, ˙x, ˙y, ˙z) = [x − x
0
, y − y
0
, z − z
0
, ˙x, ˙y, ˙z]
p
1
0 0 0 0 0
0 p
2
0 0 0 0
0 0 p
3
0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
x − x
0
y − y
0
z − z
0
˙x
˙y
˙z
(1.6)
其中：(x, y, z) 表示飞行器在 oxyz 下的位置，(x
0
, y
0
, z
0
) 表示平移末端点，(x
0
, y
0
, 0) 为安全着
陆点，z
0
为指定高度 (z
0
= 100 − 30)，( ˙x, ˙y, ˙z) = (v
x
, v
y
, v
z
) 表示速度 v 的方向分量。在平移末
端期望该状态量为 0。飞行器的 (x, y, z, ˙x, ˙y, ˙z) 由系统状态给出，x
0
, y
0
由安检系统给出，z
0
事
先给定，p
1
, p
2
, p
3
∈ R
+
，函数 ϕ
p
越大，则距离着陆点 (x
0
, y
0
) 越远，速度偏差越大。
势危函数为
ϕ
s
(x, y) =
n
i=1
|z
i
|e
−
(x−x
i
)
2
+(y−y
i
)
2
σ
2
(1.7)
其中：x, y 表示娥 3 的 x, y 轴位置，(x
i
, y
i
, z
i
) 为障碍物位置、高度/深度，i 为障碍物需要，共
n 个障碍物，σ 为超参，代表了危险威胁范围，σ 越大，障碍物对周围的危险越大。ϕ
s
越大，则
越危险。
上面给出了能量函数 (1.6) 和势危函数 (1.7)，由此，我们可以设置李雅普诺夫函数
v = ϕ(x, y, z, ˙x, ˙y, ˙z) = ϕ
p
+ k
1
ϕ
s
其中：k
1
为权重参数。由 ϕ
p
, ϕ
s
为正定，有 ϕ > 0，并且 ϕ 为状态 x(t) 的函数。由李雅普诺夫
第二定理可知：当 ϕ < 0 时，则系统在状态空间平衡点出的平衡状态是一致渐进稳定的，该平衡
点记为 v 的极小点。
下面给出李雅普诺夫最优控制律。李雅普诺夫函数的导数为三个方向周速度的函数，并且负
定。此方法可以保证娥 3 到达平移末端时，其水平速度为 0，即令
˙
ϕ = −k
x
˙x
2
− k
y
˙y
2
− k
z
˙z
2
，其
中：k
x
, k
y
, k
z
为正数。由于 ϕ 为 (x, y, z ˙x, ˙y, ˙z) 的函数，所以
dϕ
dt
可以由链式法则求出
dϕ
dt
=
∂ϕ
∂x
∂x
∂t
+
∂ϕ
∂y
∂y
∂t
+
∂ϕ
∂z
∂z
∂t
+
∂ϕ
∂ ˙x
∂ ˙x
∂t
+
∂ϕ
∂ ˙y
∂ ˙y
∂t
+
∂ϕ
∂ ˙z
∂ ˙z
∂t
=
∂ϕ
∂x
˙x +
∂ϕ
∂y
˙y +
∂ϕ
∂z
˙z +
∂ϕ
∂ ˙x
¨x +
∂ϕ
∂ ˙y
¨y +
∂ϕ
∂ ˙z
¨z
http://www.ma-xy.com 45 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
下面，我们来求解
∂ϕ
∂x
,
∂ϕ
∂ ˙x
,
∂ϕ
∂y
,
∂ϕ
∂ ˙y
,
∂ϕ
∂z
,
∂ϕ
∂ ˙z
∂ϕ
∂x
= 2p
1
(x − x
0
) −
2
σ
2
n
i=1
k
1
|z
i
|(x − x
i
)e
−
(x−x
i
)
2
+(y−y
i
)
2
σ
2
∂ϕ
∂y
= 2p
2
(y − y
0
) −
2
σ
2
n
i=1
k
1
|z
i
|(y − y
i
)e
−
(x−x
i
)
2
+(y−y
i
)
2
σ
2
∂ϕ
∂z
= 2p
3
(z − z
0
)
∂ϕ
∂ ˙x
= 2 ˙x
∂ϕ
∂ ˙y
= 2 ˙y
∂ϕ
∂ ˙z
= 2 ˙z
注：含月球自转的着陆器质心运动方程为
¨r = u + µ − 2w ˙r − w(w(r + ρ))
其中：r 为位置矢量，u 为旋转角速度，u 为控制量，µ = G
0
M。月球自转角速度 w = 2 .66 ×
10
−6
rad/s 所带来的向心力和科氏力较小，所以我们前面将其忽略了。由于着陆末端飞跨的纬度
较小，所以设月球引力场为均匀引力场。g =
G
0
M
r
2
为常量。
下面，我们来求 ( ˙x, ˙y, ˙z, ¨x, ¨y, ¨z)。( ˙x, ˙y, ˙z, ¨x, ¨y, ¨z) 即为质心运动方程，由前面的 3D 直角坐标
系下的动力系统给出。前面的 3D 直角坐标系下的动力系统为
˙x = v
x
˙y = v
y
˙z = v
z
¨x = ˙v
x
=
F
1
cos ϕ
m
−
µ
r
2
+
v
2
y
+ v
2
z
r
+
F
2x
m
¨y = ˙v
y
=
F
1
sin ϕ cos φ
m
−
v
x
v
y
+
v
2
z
r tan β
+
F
2y
m
¨z = ˙v
z
=
F
1
sin ϕ sin φ
m
−
v
x
v
y
−
v
y
v
z
r tan β
+
F
2z
m
˙m = −
F
1
+ F
2
m
其中：F
2x
, F
2y
, F
2z
为姿态发动机推力 F
2
在个方向上的分量。由此，我们的得到
dϕ
dt
。下面我们
对运动方程进行简化。
由于在着陆末端 [100 − 0]m 处，β →
π
2
, tan β → ∞，并且，要考虑姿态发动机的调整作用，
我们对娥 3 重新进行受力分析，以简化 ¨x, ¨y, ¨z。
娥 3 在末端主要有 3 个作用力：¬主发动机推力 F
1
；姿态发动机调整力 F
2
；月球引力
G
0
Mm
r
2
= µm。我们将 3 者的合力设为 F = F
1
+ F
2
+ G，设 F 在月心坐标系下各方向轴所产生
http://www.ma-xy.com 46 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
的加速度为 a
x
, a
y
, a
z
，则
¨x = a
x
= F
x
/m = (
F
1x
+
F
2x
+
G
x
)/m
¨y = a
y
= F
y
/m = (
F
1y
+
F
2y
+
G
y
)/m
¨z = a
z
= F
z
/m = (
F
1z
+
F
2z
+
G
z
)/m
为使 ˙x = −k
x
˙x
2
− k
y
˙y
2
+ k
z
˙z
2
恒成立，可以得到 x, y, z 各方向轴的期望加速度 ¨x, ¨y, ¨z 为
¨x = −k
x
˙x
2
+
n
i=1
k
1
|z
i
|e
−
(x−x
i
)
1
+(y−y
i
)
2
σ
2
x − x
i
σ
2
− p
1
(x − x
0
)
¨y = −k
y
˙x
2
+
n
i=1
k
1
|z
i
|e
−
(−x
i
)
1
+(y−y
i
)
2
σ
2
y − y
i
σ
2
− p
2
(y − y
0
)
¨z =
µ
r
− k
z
˙z
2
− p
3
(z − z
0
)
其中：k
1
, k
x
, k
y
, k
z
, p
1
, p
2
, p
3
为参数。下面我们来简化运动方程：如果按上面的运动方程进行求
解，我们还需要得到控制量 u = [F, ψ, φ, F
2x
, F
2y
, F
2z
]
T
的李雅普诺夫最优控制律。这样就太麻
烦了。
考虑到末端 β →
π
2
，并且飞行器娥 3 并不会偏离 oy 轴太远，即 α →
π
2
，月球引力常量 G
可进行化简
G ≈ [G
x
, G
y
, G
z
]
T
= [G cos β, G sin β cos α, G sin β sin α]
T
有
lim
α→π/2
β→π/2
G(α, β) = [G
x
, G
y
, G
z
]
T
= [0, 0, G]
T
说明：上式说明月球引力常量 G 只牵引娥 3 沿 z 轴向月心运动。
并且，我们考虑到主发动机推理 F
1
与 z 轴方向重合 (ψ = π/2, φ = π/2)，且姿态发动机在
娥 3 周围 (底部和上部没有)，于是
F
1x
= 0, F
1y
= 0, F
1z
= F
1
F
2z
= 0
说明：上式说明主发动机 F
1
方向与 z 轴重合，即发动机主要指引娥 3 进行 z 方向的降落；姿态
发动机 F
2
在 z 方向合力为 0，主要指引娥 3 进行 xoy 平移。所以，简化后的运动方程为
¨x = a
x
= F
2x
/m
¨y = a
y
= F
2y
/m
¨
z
=
a
z
=
F
1z
/
m
−
G
/
m
=
F
1z
/
m
−
µ
因此，化简后的控制律记为 u = [F
1
, F
2x
, F
2y
]
T
，可以由期望控制律 (期望加速度)¨x, ¨y, ¨z 求出。该
控制律 u 是在李雅普诺夫控制理论下求取的，故其可以保证李雅普诺夫稳定。
http://www.ma-xy.com 47 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
飞行器娥 3 在月面 100m 处悬停，并进行拍照，x
0
= (0, 0, 100), v
0
= (0, 0, 0)。取能量函数
加权系数 p
1
= p
2
= p
3
= 0.002，势危函数权重 k
1
= 0.1，形状系数 σ
2
= 2，势函数下降函数
k
x
= k
y
= k
z
= 0.2。
todo: 结果图：能量函数图、势危函数图、李氏函数图、在李高程图中避障轨迹、v
x
, v
y
, v
z
, x, y, z
曲线。
讨论：z
i
的确定；注：PUPF 调节器、PID 神经网络控制 (假设采用连续姿态推理)，当然，
如果不考虑姿态发动机，以 [F, ψ, φ]
T
替代亦可。
(5) 缓速下降阶段
缓速下降阶段的主要任务是控制飞行器在 4m 时的速度为 0，之后关闭发动机，娥 3 以自由
落体着陆到指定点，其示意图如图 (1.22) 所示
图 1.22: 缓速下降阶段示意图
如果娥 3 在 30m 出时，能准确降落到 P 点 (P 点位着陆点 o 的正上方 30m 处)，那么缓速
下降的模型可以简化为古老的 (经典) 软着陆问题。例如：《现代控制理论及应用》齐晓慧 P155。
缓速下降阶段的运动方程为
˙z = v
˙v =
µ
r
2
−
F
m
˙m = −
F
v
e
其中：z 为月心坐标系下的 z 轴值，即高度值，v 为速度。并且假设姿态调整力 F
2
为 0，推力
F
1
方向与 z 轴重合。
上述系统的状态量为 x = [z, v, m]
T
，控制量 u = F 。建立以燃料消耗最小为性能指标的积
分型最优控制模型如下
min
u
J =
t
f
t
0
˙m(t)dt
s.t.
x(0) = x
0
˙x = f (x, u, t)
x(t
f
) = x
f
u
min
⩽ u ⩽ u
max
http://www.ma-xy.com 48 http://www.ma-xy.com
http://www.ma-xy.com
第一章 嫦娥 3 号 1.2 嫦娥 3 号软着陆的优化控制
但是，由于实际情况 (误差原因) 飞行器在 30m 时不能达到 P 点而是达到 P
′
，因此，我们
仍要在 3 维空间中进行分析，并且还要考虑姿态发动机的存在，以实现高精度着陆工作。
承继前面的 F
1
和 F
2
的假设，假设 F
1
的方向与 z 轴重合，F
2
仅产生 F
2x
, F
2y
两个方向上
的分力，以简化模型。我们有简化后的系统动力学模型
˙x = v
x
˙y = v
y
˙z = v
z
˙v
x
= F
x
/m
˙v
y
= F
y
/m
˙v
z
= µ/r
2
− F
1
/m
˙m =
F
1
+ F
2
v
e
最优控制目标仍然为燃料消耗最少，可用适当的方法进行求解。
注：哈工大. 张锋.2009 硕士论文中给出了姿态联合动力系统，并用参数化 + 强化技术 + 约
束变换进行最优控制律的求解，并于着陆末端，在姿态耦合系统的基础上求解系统的联合镇定控
制律。
1.2.6 问题三的分析与求解
问题的分析
对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
http://www.ma-xy.com 49 http://www.ma-xy.com
http://www.ma-xy.com
1.2 嫦娥 3 号软着陆的优化控制 第一章 嫦娥 3 号
http://www.ma-xy.com 50 http://www.ma-xy.com
http://www.ma-xy.com
参考文献
[1] 陈传淼. 科学计算概率.
[2] 胡健伟, 冯怀民. 微分方程数值方法.
[3] 李乃成, 梅立泉. 数值分析.
[4] 孙军伟. 基于 SQP 方法的常推力月球软着陆轨道优化方法. 哈工大. 宇航学报.
[5] 赵吉松. 月球轨道优化快速时间逼近.
[6] 周净扬. 月球探测器软着陆精确建模及最优轨道设计. 哈工大. 宇航学报.
51