跳转至

变分法与 Euler-Lagrange 方程

前置自测

答不出 2 题以上,请先回顾线性代数与实分析基础。

  1. 什么是函数空间中的范数?\(C^1[a,b]\) 空间上的范数如何定义?
  2. 分部积分公式 \(\int_a^b u\,dv = uv\big|_a^b - \int_a^b v\,du\) 中,如果 \(u(a)=u(b)=0\),边界项会怎样?
  3. 给定 \(f:\mathbb{R}^n \to \mathbb{R}\),什么条件下 \(x^*\) 是极小值点?一阶必要条件和二阶条件分别是什么?
  4. Lagrange 乘子法的几何直觉是什么?为什么约束优化的最优点上梯度与约束面法向平行?
  5. 什么是边值问题(BVP)?它与初值问题(IVP)有何本质区别?

本章目标

学完本章后,你将能够:从第一变分为零推导任意 Lagrangian 的 Euler-Lagrange 方程并理解每一步的数学与物理意义;掌握等周约束、holonomic 约束等约束变分问题的处理方法;利用 Beltrami 恒等式与 Noether 定理快速提取守恒量简化问题求解;完整求解最速降线、悬链线等经典变分问题;理解 Lagrange 力学、测地线、minimum-snap 轨迹的变分本质;建立从变分法到 Pontryagin 极大值原理的逻辑桥梁。

知识树

变分法与 Euler-Lagrange 方程
├── 1. 泛函与变分的基本概念 ⭐
│   ├── 泛函的定义与函数空间
│   ├── 弱变分与强变分
│   └── Gateaux 导数与 Frechet 导数
├── 2. Euler-Lagrange 方程的完整推导 ⭐
│   ├── 固定端点变分问题
│   ├── 第一变分的计算
│   ├── 变分引理(du Bois-Reymond)
│   └── EL 方程的多维推广
├── 3. 经典案例完整求解 ⭐⭐
│   ├── 最速降线(Brachistochrone)
│   ├── 悬链线(Catenary)
│   ├── 测地线方程
│   └── 最小旋转面
├── 4. 带约束的变分问题 ⭐⭐
│   ├── 等周问题(Isoperimetric)
│   ├── Holonomic 约束
│   ├── 非完整约束
│   └── 横截条件(自由端点)
├── 5. Hamilton 原理与经典力学 ⭐⭐
│   ├── 最小作用量原理的真相
│   ├── Lagrange 力学与广义坐标
│   ├── Legendre 变换与 Hamilton 力学
│   └── 相空间与辛结构
├── 6. 第一积分与守恒量 ⭐⭐
│   ├── Beltrami 恒等式
│   ├── 能量守恒与动量守恒
│   └── Noether 定理
├── 7. 二阶变分与充分条件 ⭐⭐⭐
│   ├── 第二变分与 Legendre 必要条件
│   ├── Jacobi 方程与共轭点
│   ├── Weierstrass E-函数
│   └── 充分条件汇总与应用
├── 8. 直接法与数值变分 ⭐⭐⭐
│   ├── Ritz 法
│   ├── Galerkin 法
│   └── 有限元联系
└── 9. 工程应用与通向最优控制的桥梁 ⭐⭐⭐
    ├── 间接法 vs 直接法
    ├── 机器人动力学的变分基础
    ├── Minimum-snap 轨迹优化
    └── EL → PMP 的推广路径

预计阅读时间

阅读方式 时间 适合谁
精读(含手推经典问题) 10-12 小时 需要完整掌握 EL 推导、充分条件理论和 Noether 定理的读者
速读(跳过二阶条件细节) 4-5 小时 已有微积分基础、需要快速建立变分法框架的读者
速查(只看经典问题和公式表) 30 分钟 遇到具体变分问题时回来查阅 EL 方程或 Beltrami 恒等式

1.1 泛函与变分的基本概念 ⭐

动机

在有限维优化中,我们求的是使 \(f(x_1, x_2, \ldots, x_n)\) 取极值的点 \(x^* \in \mathbb{R}^n\)。但在物理和工程中,大量问题要求找的不是一个"点",而是一条"曲线"或一个"函数"。例如:两点之间哪条路径使小球下滑时间最短?哪条曲线旋转后得到的曲面面积最小?机器人关节从 A 到 B 的哪条轨迹使消耗能量最少?

这类问题的共同特征是:优化的对象不再是有限个实数,而是一个函数 \(y(x)\) 整体。我们需要一套"对函数求极值"的微积分——这就是变分法。

如果没有变分法会怎样

假设你要找两点间使小球最快滑落的曲线。如果没有变分法,你只能做以下事情:猜一条曲线(比如直线、抛物线、圆弧),计算滑落时间,然后再猜另一条,比较哪个更快。这等价于在无穷维空间中做"随机搜索"——当搜索空间是所有连续函数时,这种方法永远无法保证找到最优解。变分法提供了一个系统的方法:把无穷维问题转化为求解微分方程,一步到位得到最优函数。

历史:从 Bernoulli 兄弟的挑战开始

1696 年 6 月,Johann Bernoulli 在《学术纪事》(Acta Eruditorum) 上发布了一道挑战题:

"给定竖直平面内两个不在同一竖直线上的点 A 和 B,求连接 A、B 的曲线,使得仅受重力作用的质点从 A 滑到 B 所需时间最短。"

这就是著名的**最速降线问题**(Brachistochrone,源自希腊语 brachistos = 最短 + chronos = 时间)。六位数学家给出了解答:Johann Bernoulli、Jakob Bernoulli、Newton、Leibniz、L'Hopital 和 Tschirnhaus。Jakob Bernoulli 的解法尤其重要——他没有使用弟弟 Johann 的光学类比技巧,而是发展了一种通用方法,通过"令变分为零"得到曲线必须满足的微分方程。

Euler 在 1744 年的《求具有某种极值性质的曲线的方法》(Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes) 中系统化了这一方法,首次写下了我们今天称为 Euler-Lagrange 方程的公式。Lagrange 在 1755 年给 Euler 的信中引入了 \(\delta\) 记号(变分符号),使推导更加优雅——Euler 读后十分赞赏,等了将近 30 年才发表自己的方法,以让年轻的 Lagrange 获得优先权。

这段历史揭示了变分法诞生的内在逻辑:

年代 人物 贡献 方法论进步
1696 Johann Bernoulli 提出最速降线问题 将物理问题数学化
1697 Jakob Bernoulli 通用变分方法 从特殊到一般的方法论飞跃
1744 Euler 系统化为 EL 方程 建立完整理论框架
1755 Lagrange 引入 \(\delta\) 记号 使计算规范化
1834 Hamilton 最小作用量原理 统一力学形式
1879 Weierstrass 充分条件理论 完善数学严格性

泛函的严格定义

定义 1.1.1(泛函):设 \(\mathcal{V}\) 是一个函数空间(如 \(C^1[a,b]\)),映射 \(J: \mathcal{V} \to \mathbb{R}\) 称为**泛函**(functional),即将一个函数映射为一个实数。

最常见的形式是积分型泛函:

\[J[y] = \int_a^b L(x, y(x), y'(x))\,dx\]

其中 \(L: \mathbb{R} \times \mathbb{R} \times \mathbb{R} \to \mathbb{R}\) 称为 Lagrangian(被积函数),\(y \in C^1[a,b]\) 是待优化的函数。

本质洞察:泛函与函数的关系,类似于函数与数的关系。函数将数映射为数 \(f: \mathbb{R} \to \mathbb{R}\);泛函将函数映射为数 \(J: \mathcal{V} \to \mathbb{R}\)。有限维优化的工具(梯度、极值条件)需要被"抬升"到无穷维空间中——这就是变分法的核心任务。

类比:如果把函数 \(f(x_1, \ldots, x_n)\) 想象成一座有限维的山,那么泛函 \(J[y]\) 就是一座无穷维的山——每个"点"是一条曲线,"海拔"是积分值。变分法就是在这座无穷维的山上找极值点的工具。

为了使类比更精确,我们列出有限维到无穷维的完整对应关系:

有限维(\(\mathbb{R}^n\) 无穷维(函数空间)
向量 \(x = (x_1, \ldots, x_n)\) 函数 \(y(x)\)
函数 \(f(x): \mathbb{R}^n \to \mathbb{R}\) 泛函 \(J[y]: \mathcal{V} \to \mathbb{R}\)
梯度 \(\nabla f\) 变分导数 \(\frac{\delta J}{\delta y}\)
\(\nabla f = 0\)(驻点条件) EL 方程(\(\delta J = 0\)
Hessian \(\nabla^2 f \succeq 0\) 第二变分 \(\delta^2 J \geq 0\)
Taylor 展开 \(f(x+h) = f(x) + \nabla f \cdot h + \frac{1}{2}h^\top H h + \cdots\) \(J[y+\varepsilon h] = J[y] + \varepsilon\delta J + \frac{\varepsilon^2}{2}\delta^2 J + \cdots\)

这个对应表是理解变分法的指南针——每当你在无穷维中迷失方向时,回到有限维找到对应概念,就能找到前进的方向。

容许函数空间

变分问题不是在所有函数中搜索,而是在满足特定条件的函数集合中搜索。

定义 1.1.2(容许函数类):对于固定端点问题 \(y(a) = A\), \(y(b) = B\),容许函数空间为

\[\mathcal{A} = \{y \in C^k[a,b] : y(a) = A,\; y(b) = B\}\]

其中 \(k\) 取决于 Lagrangian 中出现的最高阶导数。对标准 \(L(x,y,y')\) 问题,取 \(k=1\)\(k=2\)

为什么要指定光滑性? 如果不限制函数的光滑程度,可能出现"锯齿形"函数,使得积分定义本身出问题(\(y'\) 可能处处不存在)。不同的光滑性假设对应不同的变分问题类:\(C^2\) 对应经典变分(可以做分部积分),\(C^1\) 允许"角点"(需要 Weierstrass-Erdmann 角条件),\(W^{1,2}\)(Sobolev 空间)是现代 PDE 变分方法的标准设定。

弱变分与强变分

函数空间上可以定义不同的范数,导致"邻域"的概念不同,进而产生不同强度的极值概念。

概念 范数 直觉 含义
弱变分 \(C^1\) 范数 $|h|_1 = \max h + \max
强变分 \(C^0\) 范数 $|h|_0 = \max h $

为什么这个区分重要? 弱极小比强极小更容易达到——一条曲线可能在所有"光滑的邻近曲线"中最优(弱极小),但在允许"急转弯"的竞争者面前就不再最优(非强极小)。EL 方程是两种极值的共同必要条件,但充分条件不同:弱极小只需 Legendre 条件 \(L_{y'y'} \geq 0\);强极小还需要 Weierstrass E-函数非负。

反事实推理:如果我们不区分弱极小和强极小会怎样?想象一条最短路径:在光滑曲线中它是最短的(弱极小),但一条来回折叠的曲线(\(C^0\) 范数小但导数大)可能在投影距离意义上更短。不区分这两种极值,就无法正确判断解的稳定性——在工程中,这意味着你的轨迹可能被"小扰动"(光滑)保持最优,却被"急转弯"(非光滑)打败。最优控制中的 bang-bang 控制正是"强变分"邻域中的竞争者。

Gateaux 导数与 Frechet 导数

定义 1.1.3(Gateaux 导数 / 第一变分):泛函 \(J\)\(y\) 处沿方向 \(h\) 的 Gateaux 导数定义为

\[\delta J[y; h] = \lim_{\varepsilon \to 0} \frac{J[y + \varepsilon h] - J[y]}{\varepsilon} = \frac{d}{d\varepsilon}\bigg|_{\varepsilon=0} J[y + \varepsilon h]\]

这就是**第一变分**。它完全类比于有限维中的方向导数 \(Df(x; v) = \lim_{\varepsilon \to 0} \frac{f(x+\varepsilon v)-f(x)}{\varepsilon}\)

定义 1.1.4(Frechet 导数):若存在连续线性泛函 \(DJ_y: \mathcal{V} \to \mathbb{R}\) 使得

\[J[y+h] = J[y] + DJ_y(h) + o(\|h\|)\]

则称 \(J\)\(y\) 处 Frechet 可微。Frechet 可微蕴含 Gateaux 可微,反之不然。

反事实推理:如果我们只使用 Frechet 导数会怎样?在很多变分问题中(特别是非线性问题),Frechet 导数可能不存在,但 Gateaux 导数(方向导数)仍然存在。变分法的历史发展表明,先计算方向导数并令其为零(\(\delta J = 0\)),然后用变分引理提取结论——这条路径不需要全导数的存在性,因此适用范围更广。

两者的关系与选择

  • Gateaux 导数:只要求沿每个方向的极限存在,不要求关于方向的一致性。它是"逐方向"检查的工具,计算简单,但一般不保证连续依赖于方向。
  • Frechet 导数:要求余项 \(o(\|h\|)\) 对所有方向一致成立。它是更强的条件,保证了全局的线性近似质量。
  • 变分法的实践:几乎所有推导只使用 Gateaux 导数(第一变分),因为我们只需要"对所有方向第一变分为零"这一条件。Frechet 导数在理论分析(如泛函的 Taylor 展开、非线性泛函分析)中更重要。

极值的必要条件

定理 1.1.1(极值的一阶必要条件):若 \(y^*\) 是泛函 \(J\) 的(弱/强)极小值,则对所有容许方向 \(h\)(满足端点条件 \(h(a)=h(b)=0\)),第一变分为零:

\[\delta J[y^*; h] = 0, \quad \forall h\]

证明:定义 \(\Phi(\varepsilon) = J[y^* + \varepsilon h]\)。若 \(y^*\) 是极小,则 \(\Phi(0) \leq \Phi(\varepsilon)\) 对足够小的 \(|\varepsilon|\) 成立。因此 \(\Phi'(0) = 0\),而 \(\Phi'(0) = \delta J[y^*; h]\)\(\square\)

这完全类比于有限维中 \(f\) 在极值点 \(x^*\) 处梯度为零 \(\nabla f(x^*) = 0\)

为什么这个类比如此重要? 它说明变分法不是一门全新的学科,而是微积分从有限维到无穷维的自然推广。你在大一学的"极值点处导数为零"这一朴素思想,经过适当的推广,就能处理函数空间中的优化问题——这体现了数学思想的统一性和延伸性。

变分导数与泛函的 Taylor 展开

有限维中,函数 \(f:\mathbb{R}^n \to \mathbb{R}\) 在点 \(x\) 处的 Taylor 展开为

\[f(x+h) = f(x) + \nabla f(x)^\top h + \frac{1}{2}h^\top \nabla^2 f(x) h + \cdots\]

泛函完全类似。对 \(J[y] = \int_a^b L(x,y,y')dx\),在 \(y^*\) 处沿方向 \(h\) 的 Taylor 展开为

\[J[y^* + \varepsilon h] = J[y^*] + \varepsilon\,\delta J[y^*; h] + \frac{\varepsilon^2}{2}\,\delta^2 J[y^*; h] + \cdots\]

这就是泛函的展开。其中第一变分 \(\delta J\) 类比于梯度 \(\nabla f\),第二变分 \(\delta^2 J\) 类比于 Hessian \(\nabla^2 f\)

如果第一变分可以写成

\[\delta J[y^*; h] = \int_a^b \frac{\delta J}{\delta y}(x) \cdot h(x)\,dx\]

那么 \(\frac{\delta J}{\delta y}(x)\) 称为**变分导数**(functional derivative 或 variational derivative)。它是"泛函的梯度"——在函数空间中指向 \(J\) 增长最快的方向。EL 方程 \(\frac{\delta J}{\delta y} = L_y - \frac{d}{dx}L_{y'} = 0\) 就是"变分导数为零"条件。

变分导数的物理意义\(\frac{\delta J}{\delta y}(x_0)\) 衡量的是"在点 \(x_0\) 处微小改变函数值 \(y(x_0)\) 时,泛函 \(J\) 的变化率"。它是逐点定义的——这和有限维中梯度的每个分量 \(\frac{\partial f}{\partial x_i}\) 衡量"改变第 \(i\) 个坐标的效果"完全类比。

练习

  1. 概念题:解释为什么"泛函的极值"不能简单地通过"对 \(y(x)\) 的每个点分别求极值"来得到。提示:考虑 \(J[y] = \int_0^1 (y')^2 dx\)\(y(0)=0, y(1)=1\)
  2. 计算题:对泛函 \(J[y] = \int_0^1 [(y')^2 + y^2]\,dx\),计算在 \(y_0(x) = x\) 处沿方向 \(h(x) = \sin(\pi x)\) 的第一变分。
  3. 思考题:在什么意义下,有限维 \(f:\mathbb{R}^n \to \mathbb{R}\) 的梯度 \(\nabla f\) 对应于无穷维中的什么概念?(提示:Riesz 表示定理)
  4. 计算题:对泛函 \(J[y] = \int_0^1 y\sqrt{1+y'^2}\,dx\)(悬链线问题),计算变分导数 \(\frac{\delta J}{\delta y}\)

1.2 Euler-Lagrange 方程的完整推导 ⭐

动机

上一节建立了极值的一阶必要条件:\(\delta J = 0\) 对所有容许变分 \(h\) 成立。但这是一个关于"所有方向 \(h\)"的条件——我们如何从中提取出 \(y^*\) 本身必须满足的微分方程?答案是:先计算第一变分的显式表达式,然后利用一个关键引理(变分引理)将积分条件转化为逐点条件。

这一步的逻辑结构是:

\[\text{对所有 } h, \int_a^b (\cdots) h\,dx = 0 \quad \xrightarrow{\text{变分引理}} \quad (\cdots) = 0, \; \forall x\]

从"积分为零"到"被积函数为零"——这个跳跃不是显然的(一个非零函数的积分可以为零,比如 \(\sin x\)\([0, 2\pi]\) 上的积分)。变分引理的深刻之处在于:它利用了"对所有 \(h\) 都成立"这一条件的全部力量——如果被积函数在某点不为零,我们总能构造一个集中在该点附近的 \(h\) 使积分不为零,从而产生矛盾。

完整推导

考虑标准问题:在 \(y(a)=A\), \(y(b)=B\) 下极小化

\[J[y] = \int_a^b L(x, y, y')\,dx\]

Step 1:参数化扰动

\(y^*\) 是极值函数。取任意 \(\eta \in C^2[a,b]\) 满足 \(\eta(a) = \eta(b) = 0\)(端点条件保证扰动后的函数仍在容许类中)。定义单参数族

\[y_\varepsilon(x) = y^*(x) + \varepsilon \eta(x)\]

\(y_\varepsilon \in \mathcal{A}\) 对所有 \(\varepsilon\)(因为端点条件被 \(\eta\) 的零边界保持)。

为什么取 \(\eta(a)=\eta(b)=0\) 因为我们考虑的是固定端点问题——\(y(a)=A\)\(y(b)=B\) 是硬约束。任何扰动都不能违反这个约束,所以扰动函数 \(\eta\) 在端点处必须为零。如果端点自由,\(\eta\) 在端点的值就不受限制——这将导致"自然边界条件",我们在 1.8 节讨论。

Step 2:计算第一变分

\[\delta J = \frac{d}{d\varepsilon}\bigg|_{\varepsilon=0} J[y^* + \varepsilon\eta]\]
\[= \frac{d}{d\varepsilon}\bigg|_{\varepsilon=0} \int_a^b L(x, y^* + \varepsilon\eta, y^{*\prime} + \varepsilon\eta')\,dx\]

由 Leibniz 积分法则(对参数 \(\varepsilon\) 求导可以和积分交换,这需要 \(L\) 的连续可微性):

\[= \int_a^b \left[\frac{\partial L}{\partial y}\bigg|_{y^*} \cdot \eta + \frac{\partial L}{\partial y'}\bigg|_{y^*} \cdot \eta'\right] dx\]

\(L_y = \frac{\partial L}{\partial y}\), \(L_{y'} = \frac{\partial L}{\partial y'}\)(在 \(y^*\) 处取值),则

\[\delta J = \int_a^b \left[L_y \cdot \eta + L_{y'} \cdot \eta'\right] dx\]

这一步的数学意义:我们把"泛函在 \(y^*\) 处的方向导数"写成了两项之和——第一项 \(L_y \eta\) 度量的是 \(L\) 对"函数值变化"的敏感度,第二项 \(L_{y'}\eta'\) 度量的是 \(L\) 对"斜率变化"的敏感度。

Step 3:分部积分

第二项 \(\int_a^b L_{y'} \cdot \eta'\,dx\) 含有 \(\eta'\),而我们希望得到只含 \(\eta\) 的表达式(因为变分引理需要 \(\int f \cdot \eta = 0\) 的形式)。对第二项做分部积分:

\[\int_a^b L_{y'} \cdot \eta'\,dx = \left[L_{y'} \cdot \eta\right]_a^b - \int_a^b \frac{d}{dx}L_{y'} \cdot \eta\,dx\]

关键步骤:由于 \(\eta(a) = \eta(b) = 0\),边界项

\[\left[L_{y'} \cdot \eta\right]_a^b = L_{y'}(b) \cdot \eta(b) - L_{y'}(a) \cdot \eta(a) = 0\]

完全消失!这就是为什么固定端点条件如此重要——它消除了边界项,使得推导可以继续。

Step 4:合并

\[\delta J = \int_a^b \left[L_y - \frac{d}{dx}L_{y'}\right] \eta\,dx = 0\]

这对所有满足 \(\eta(a)=\eta(b)=0\)\(\eta \in C^2[a,b]\)(甚至 \(C_0^\infty(a,b)\))成立。

Step 5:变分引理(du Bois-Reymond 引理)

这是从积分条件到逐点条件的关键桥梁。

引理 1.2.1(变分基本引理 / du Bois-Reymond 引理):若 \(f \in C[a,b]\),且对所有 \(\eta \in C_0^\infty(a,b)\)

\[\int_a^b f(x) \eta(x)\,dx = 0\]

\(f(x) \equiv 0\)\([a,b]\) 上恒成立。

证明:反设 \(f(x_0) > 0\) 对某个 \(x_0 \in (a,b)\)\(f < 0\) 类似)。由连续性,存在 \(\delta > 0\) 使得 \(f(x) > 0\)\((x_0-\delta, x_0+\delta) \subset (a,b)\) 上。构造

\[\eta(x) = \begin{cases} \exp\left(-\frac{1}{(\delta^2 - (x-x_0)^2)}\right), & |x-x_0| < \delta \\ 0, & \text{otherwise} \end{cases}\]

这是一个 bump 函数:\(\eta \in C_0^\infty(a,b)\)\(\eta \geq 0\),且在 \((x_0-\delta, x_0+\delta)\) 上严格正。于是

\[\int_a^b f\eta\,dx = \int_{x_0-\delta}^{x_0+\delta} f(x)\eta(x)\,dx > 0\]

矛盾。故 \(f \equiv 0\)\(\square\)

为什么需要这个引理的证明如此仔细? 因为变分引理是整个变分法的逻辑基石。没有它,我们无法从"积分为零"跳到"被积函数为零"——而这个跳跃正是把无穷维条件(\(\delta J = 0\) 对所有 \(h\))转化为有限维条件(逐点 ODE)的关键。数学上,这个引理的本质是 \(C_0^\infty\)\(L^2\) 中的稠密性——测试函数族足够丰富,可以"探测"被积函数在任何一点的非零性。

变分引理的更强版本(du Bois-Reymond):上述引理要求 \(f\) 连续。du Bois-Reymond 1879 年证明了更强的版本:即使只假设 \(f\) 可积,\(\int f\eta = 0\) 对所有 \(\eta \in C_0^\infty\) 成立也蕴含 \(f\) 几乎处处为零。这个加强版在现代 PDE 弱解理论中至关重要。

Step 6:得到 Euler-Lagrange 方程

由变分引理,从 \(\int_a^b [L_y - \frac{d}{dx}L_{y'}]\eta\,dx = 0\) 对所有 \(\eta\) 成立,得

\[\boxed{\frac{\partial L}{\partial y} - \frac{d}{dx}\frac{\partial L}{\partial y'} = 0}\]

这就是**Euler-Lagrange 方程**——变分法最核心的结果。

EL 方程的物理意义

EL 方程可以改写为

\[\frac{d}{dx}\frac{\partial L}{\partial y'} = \frac{\partial L}{\partial y}\]

左边 \(\frac{d}{dx}L_{y'}\) 是"广义动量的变化率",右边 \(L_y\) 是"广义力"。这完全类比于 Newton 第二定律 \(\dot{p} = F\)。事实上,对力学系统取 \(L = T - V\),EL 方程就是 Newton 方程的等价形式。

跨领域类比:EL 方程之于泛函,如同梯度为零之于函数。就像 \(\nabla f = 0\) 把"找极值点"转化为"解代数方程",EL 方程把"找极值函数"转化为"解微分方程"。这个降维——从无穷维搜索到解 ODE——是变分法的核心威力。

本质洞察:Euler-Lagrange 方程本质上是一个"无穷维的梯度为零条件"。就像有限维中 \(\nabla f = 0\) 给出极值点满足的代数方程组一样,\(\delta J = 0\) 通过 EL 方程给出极值函数满足的微分方程。变分法把"在函数空间中搜索"降维为"解微分方程"——这是其威力所在。

展开形式

\(\frac{d}{dx}L_{y'}\) 用链式法则展开:

\[\frac{d}{dx}L_{y'} = L_{y'x} + L_{y'y} \cdot y' + L_{y'y'} \cdot y''\]

因此 EL 方程写成

\[L_y - L_{y'x} - L_{y'y} \cdot y' - L_{y'y'} \cdot y'' = 0\]

这是关于 \(y\) 的**二阶常微分方程**(假设 \(L_{y'y'} \neq 0\))。配合两个端点条件 \(y(a)=A\), \(y(b)=B\),构成一个两点边值问题(BVP)。

如果 \(L_{y'y'} = 0\) 会怎样? 这意味着 \(L\) 关于 \(y'\) 是线性的(或常数),EL 方程退化为一阶甚至零阶方程。这种退化情形在最优控制中对应"奇异弧",是 PMP 理论中的重要课题。

EL 方程的正则性问题:我们的推导假设极值函数 \(y^* \in C^2\)。如果 \(y^*\) 只是 \(C^1\)(存在"角点",即 \(y''\) 在某些点不连续),EL 方程在角点不成立,取而代之的是 Weierstrass-Erdmann 角条件\(L_{y'}\)\(L - y'L_{y'}\) 在角点处连续。这在光学(折射定律)和最优控制(切换点)中有重要应用。

多维推广

当未知函数有多个分量 \(\mathbf{y}(x) = (y_1(x), \ldots, y_n(x))\) 时,泛函变为

\[J[\mathbf{y}] = \int_a^b L(x, y_1, \ldots, y_n, y_1', \ldots, y_n')\,dx\]

EL 方程变为 \(n\) 个耦合方程的系统:

\[\frac{\partial L}{\partial y_i} - \frac{d}{dx}\frac{\partial L}{\partial y_i'} = 0, \quad i = 1, \ldots, n\]

这正是 Lagrange 力学的标准形式:取 \(x = t\)(时间),\(\mathbf{y} = \mathbf{q}\)(广义坐标),\(L = T - V\)(动能减势能),则 EL 方程给出运动方程。

向量形式:用向量记号,设 \(\mathbf{q} \in \mathbb{R}^n\),则 EL 方程可以紧凑地写为

\[\frac{d}{dt}\nabla_{\dot{\mathbf{q}}} L - \nabla_{\mathbf{q}} L = \mathbf{0}\]

\(L = \frac{1}{2}\dot{\mathbf{q}}^\top M(\mathbf{q})\dot{\mathbf{q}} - V(\mathbf{q})\)(机器人系统的标准 Lagrangian),展开后恰好给出

\[M(\mathbf{q})\ddot{\mathbf{q}} + C(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \mathbf{0}\]

其中 Coriolis 矩阵 \(C\) 的元素是 \(M(\mathbf{q})\) 在构型空间上的 Christoffel 符号——这不是巧合,而是因为 Lagrange 力学本质上就是构型流形上的测地线方程加上势能项。

高阶情形

当 Lagrangian 含高阶导数 \(L(x, y, y', y'', \ldots, y^{(n)})\) 时,通过反复分部积分得到 Euler-Poisson 方程

\[\sum_{k=0}^{n} (-1)^k \frac{d^k}{dx^k} \frac{\partial L}{\partial y^{(k)}} = 0\]

\(n=2\)(含二阶导数的 Lagrangian):

\[\frac{\partial L}{\partial y} - \frac{d}{dx}\frac{\partial L}{\partial y'} + \frac{d^2}{dx^2}\frac{\partial L}{\partial y''} = 0\]

这是四阶 ODE,需要四个边界条件。工程应用:minimum-jerk 轨迹优化 \(L = (\dddot{x})^2\) 给出六阶 ODE,解为五次多项式;minimum-snap 优化 \(L = (x^{(4)})^2\) 给出八阶 ODE,解为七次多项式——这正是四旋翼轨迹规划的基础。

推导 Euler-Poisson 方程的关键步骤(以 \(n=2\) 为例):

第一变分为

\[\delta J = \int_a^b \left[L_y \eta + L_{y'}\eta' + L_{y''}\eta''\right]dx\]

\(L_{y'}\eta'\) 做一次分部积分,对 \(L_{y''}\eta''\) 做两次分部积分:

\[\int L_{y'}\eta'\,dx = [L_{y'}\eta]_a^b - \int (L_{y'})'\eta\,dx\]
\[\int L_{y''}\eta''\,dx = [L_{y''}\eta']_a^b - [L_{y''}'\eta]_a^b + \int (L_{y''})'' \eta\,dx\]

\(\eta(a) = \eta(b) = \eta'(a) = \eta'(b) = 0\) 时(四个边界条件),所有边界项消失,合并得

\[\delta J = \int_a^b [L_y - (L_{y'})' + (L_{y''})'']\eta\,dx = 0\]

由变分引理得到四阶 Euler-Poisson 方程。符号的交替 \((-1)^k\) 来自每次分部积分引入的一个负号。

⚠️ 常见陷阱

陷阱 1(概念误区):认为 EL 方程给出的是极小值

  • 新手想法:"EL 方程的解就是使 \(J\) 最小的函数"
  • 实际情况:EL 方程只是必要条件(类似有限维中 \(\nabla f = 0\)),其解可能是极小值、极大值或鞍点
  • 正确做法:求出 EL 方程的解后,必须用二阶条件(Legendre、Jacobi)验证是否为极小
  • 自检方法:对简单问题,可以微扰 EL 解并计算泛函值,看是增大还是减小

陷阱 2(推导错误):分部积分时忘记边界项

  • 错误做法:直接写 \(\int L_{y'}\eta' = -\int (L_{y'})'\eta\) 不检查边界
  • 根本原因:只有当 \(\eta\) 在端点为零时边界项才消失;自由端点时边界项给出自然边界条件
  • 正确做法:永远先写出完整的分部积分公式 \(\int u\,dv = [uv]_a^b - \int v\,du\),然后根据端点条件判断边界项
  • 延伸:对高阶问题,多次分部积分时边界项更复杂,必须逐步仔细跟踪

陷阱 3(技术错误):对 \(\frac{d}{dx}L_{y'}\) 的展开搞错

  • 错误做法:以为 \(\frac{d}{dx}L_{y'} = L_{y'y'}y''\),漏掉 \(L_{y'x}\)\(L_{y'y}y'\)
  • 根本原因:\(L_{y'}\)\(x\)\(y(x)\)\(y'(x)\) 三者的复合函数,链式法则给出三项
  • 正确做法:完整展开 \(\frac{d}{dx}L_{y'}(x, y(x), y'(x)) = L_{y'x} + L_{y'y}y' + L_{y'y'}y''\)

陷阱 4(概念误区):混淆 \(\frac{\partial}{\partial y}\)\(\frac{d}{dx}\)

  • 错误想法:"\(L_y\)\(L\)\(y\) 求导后再代入 \(y(x)\)"
  • 实际要求:\(L_y\) 是把 \(L(x,y,p)\) 视为三个**独立变量**的函数,对第二个变量 \(y\) 求偏导,然后才代入 \(y = y^*(x)\), \(p = y^{*\prime}(x)\)
  • 这个区分在 Hamilton 形式中更加关键:\(H(q,p,t)\)\(q\)\(p\) 是独立变量

练习

  1. \(J[y] = \int_0^1 (y'^2 + 2yy')\,dx\)\(y(0)=0\), \(y(1)=1\),写出 EL 方程并求解。
  2. 证明:如果 \(L\) 不含 \(y\)(即 \(L = L(x, y')\)),则 EL 方程退化为 \(\frac{d}{dx}L_{y'} = 0\),即 \(L_{y'} = \text{const}\)
  3. 跨章综合题:对二维问题 \(J[y_1, y_2] = \int_0^T [(\dot{y}_1)^2 + (\dot{y}_2)^2 + y_1^2 + y_2^2]\,dt\),写出耦合 EL 方程系统,并解释其物理意义(提示:这是两个耦合谐振子)。

1.3 经典案例完整求解 ⭐⭐

动机

EL 方程是一个通用的工具——但要真正掌握它,必须在具体问题上反复练习。本节精选四个在数学史和工程中最重要的变分问题,从头到尾完整求解,展示 EL 方程与各种简化技巧的配合使用。

1.3.1 最速降线问题(Brachistochrone)

问题陈述:在竖直平面内,质点从 \(A = (0, 0)\) 仅受重力沿曲线 \(y(x)\) 滑到 \(B = (x_1, y_1)\)(取 \(y\) 轴向下为正),求使滑行时间最短的曲线。

为什么这个问题如此重要? 它是人类第一次系统地研究"最优曲线"问题(1696年),其解决方案催生了整个变分法的诞生。此外,"先陡后缓"的最速降线策略在工程中有深刻意义——机器人轨迹规划中"先加速后匀速"的策略、快递无人机的最优航路规划等都体现了同样的原理。

建模:由能量守恒,速度 \(v = \sqrt{2gy}\)\(y\) 向下为正,初速度为零)。弧长微元 \(ds = \sqrt{1 + y'^2}\,dx\),时间 \(dt = ds/v\),故

\[T[y] = \int_0^{x_1} \frac{\sqrt{1 + y'^2}}{\sqrt{2gy}}\,dx\]

Lagrangian 为 \(L = \sqrt{\frac{1 + y'^2}{2gy}}\)

求解策略:观察到 \(L\) 不显含 \(x\)(即 \(\partial L / \partial x = 0\)),可以使用 **Beltrami 恒等式**直接得到一个一阶方程,而不必展开完整的二阶 EL 方程。

Beltrami 恒等式(详见 1.6 节):当 \(L = L(y, y')\) 不显含 \(x\) 时,

\[L - y' \frac{\partial L}{\partial y'} = C \quad (\text{常数})\]

计算 \(L_{y'}\)

\[\frac{\partial L}{\partial y'} = \frac{1}{\sqrt{2gy}} \cdot \frac{y'}{\sqrt{1+y'^2}}\]

代入 Beltrami 恒等式

\[\frac{\sqrt{1+y'^2}}{\sqrt{2gy}} - y' \cdot \frac{y'}{\sqrt{2gy}\sqrt{1+y'^2}} = C\]
\[\frac{1}{\sqrt{2gy}} \left[\frac{1+y'^2 - y'^2}{\sqrt{1+y'^2}}\right] = C\]
\[\frac{1}{\sqrt{2gy}\sqrt{1+y'^2}} = C\]

两边平方:

\[\frac{1}{2gy(1+y'^2)} = C^2\]

整理得:

\[y(1 + y'^2) = \frac{1}{2gC^2} =: 2r\]

其中 \(r\) 为待定常数。

参数化求解:令 \(y' = \cot(\theta/2)\)(使 \(1 + y'^2 = 1 + \cot^2(\theta/2) = \csc^2(\theta/2)\)),代入 \(y \cdot \csc^2(\theta/2) = 2r\),得 \(y = 2r\sin^2(\theta/2) = r(1 - \cos\theta)\)

为什么选 \(y' = \cot(\theta/2)\) 作为参数化? 这不是凭空而来的——观察方程 \(y(1+y'^2) = 2r\),如果令 \(y = 2r\sin^2(\theta/2)\)(受到三角恒等式 \(1+\cot^2\alpha = \csc^2\alpha\) 的启发),则 \(1+y'^2 = 2r/y = \csc^2(\theta/2)\),因此 \(y'^2 = \csc^2(\theta/2) - 1 = \cot^2(\theta/2)\),给出 \(y' = \cot(\theta/2)\)。参数化的目的是把隐式 ODE 转化为可以直接积分的形式。

\(y' = dy/dx = \cot(\theta/2)\)\(dy = r\sin\theta\,d\theta\)

\[dx = \frac{dy}{y'} = \frac{r\sin\theta}{\cot(\theta/2)}\,d\theta = r\sin\theta \cdot \frac{\sin(\theta/2)}{\cos(\theta/2)}\,d\theta\]

利用 \(\sin\theta = 2\sin(\theta/2)\cos(\theta/2)\)

\[dx = r \cdot 2\sin^2(\theta/2)\,d\theta = r(1 - \cos\theta)\,d\theta\]

积分得 \(x = r(\theta - \sin\theta) + \text{const}\)。由 \(A = (0,0)\) 对应 \(\theta = 0\),常数为零。

最终结果:最速降线是**摆线**(cycloid)

\[\boxed{x = r(\theta - \sin\theta), \quad y = r(1 - \cos\theta)}\]

这恰好是半径为 \(r\) 的圆在水平线上滚动时,圆周上一点的轨迹。

物理直觉:为什么最快的路径不是直线?因为曲线初始段陡峭→质点迅速获得高速,然后以高速运行较长的水平距离。虽然路径更长,但"先加速再匀速"的策略总时间比"匀加速走直线"更短。这个原理在机器人轨迹规划中同样适用:有时绕路但走得快,比走直线但慢更划算。

确定常数 \(r\):常数 \(r\) 由终点条件 \(B = (x_1, y_1)\) 确定。需要解

\[x_1 = r(\theta_1 - \sin\theta_1), \quad y_1 = r(1 - \cos\theta_1)\]

这是一个关于 \((r, \theta_1)\) 的超越方程组,一般无解析解,需数值求解。但对特殊情形(如 \(B\) 恰好在一条完整摆线拱上),有解析解。

数值求解策略:从第二个方程 \(y_1 = r(1-\cos\theta_1)\)\(r = y_1/(1-\cos\theta_1)\),代入第一个方程:

\[x_1 = \frac{y_1(\theta_1 - \sin\theta_1)}{1 - \cos\theta_1}\]

这是关于 \(\theta_1\) 的单变量超越方程,可用 Newton 法或二分法高效求解。

代价对比:对 \(B = (1, 1)\)\(g = 9.81\)),用数值方法可以计算: - 直线路径时间:\(T_{\text{line}} \approx 0.6389\,\text{s}\) - 摆线(最速降线)时间:\(T_{\text{cycloid}} \approx 0.5816\,\text{s}\) - 时间节省约 \(9\%\)

看起来不多——但这是全局最优!对于更极端的几何(\(B\) 远低于 \(A\) 且水平距离大),时间节省可达 \(20\%\) 以上。

1.3.2 悬链线问题(Catenary)

问题陈述:均匀密度的柔性绳索,两端固定在 \((a, y_a)\)\((b, y_b)\),在重力下的平衡形状是什么?

历史背景:Galileo (1638) 曾猜测悬挂的链条形状是抛物线。Jungius (1669) 用实验证明不是抛物线。Huygens、Leibniz、Johann Bernoulli 在 1691 年几乎同时正确求解了这个问题——答案是双曲余弦函数。

建模:绳索在重力场中的势能正比于 \(\int y\,ds = \int y\sqrt{1+y'^2}\,dx\)\(y\) 向上为正)。平衡态使势能最小,故

\[J[y] = \int_a^b y\sqrt{1+y'^2}\,dx\]

Lagrangian \(L = y\sqrt{1+y'^2}\) 不显含 \(x\)

用 Beltrami 恒等式

\[L - y'L_{y'} = y\sqrt{1+y'^2} - y' \cdot \frac{yy'}{\sqrt{1+y'^2}} = \frac{y}{\sqrt{1+y'^2}} = c\]

因此 \(y^2 = c^2(1+y'^2)\),即 \(y'^2 = (y/c)^2 - 1\)

分离变量:\(\frac{dy}{\sqrt{(y/c)^2 - 1}} = dx\)。令 \(y = c\cosh u\)\(dy = c\sinh u\,du\)

\[\frac{c\sinh u\,du}{\sqrt{\cosh^2 u - 1}} = \frac{c\sinh u\,du}{\sinh u} = c\,du = dx\]

因此 \(u = (x - x_0)/c\),故

\[\boxed{y(x) = c\cosh\frac{x - x_0}{c}}\]

这就是**悬链线**。参数 \(c\)\(x_0\) 由两端固定条件确定。

悬链线与抛物线的对比

性质 悬链线 \(y = c\cosh(x/c)\) 抛物线 \(y = ax^2 + bx + c\)
物理来源 均匀密度绳在重力下的平衡 均匀载荷(如桥面自重)下的平衡
二阶近似 在顶点附近 \(\cosh u \approx 1 + u^2/2\),与抛物线一致 精确抛物线
远离顶点 指数增长 $e^{ x/c
工程应用 无载荷绳索/电线 桥梁主缆(载荷均匀分布在水平方向)

这解释了 Galileo 的错误:在绳索下垂较小时(\(x/c \ll 1\)),\(\cosh\) 和抛物线几乎不可区分——只有在大下垂比时差异才明显。

工程应用:悬链线出现在桥梁悬索、电力线路、绳驱动机器人的绳索建模中。反向悬链线(倒挂的悬链线拱)是最优的石拱形状——受压而无弯矩,这是 Hooke 在 1675 年发现的("As hangs a flexible chain, so but inverted stands the rigid arch")。

1.3.3 测地线方程

问题陈述:在 Riemann 流形 \((M, g)\) 上,连接两点的最短曲线(测地线)满足什么方程?

建模:曲线 \(\gamma(t) = (x^1(t), \ldots, x^n(t))\) 的长度为

\[\mathcal{L}[\gamma] = \int_0^1 \sqrt{g_{ij}(\gamma)\dot{x}^i\dot{x}^j}\,dt\]

但由于长度泛函具有参数化不变性(重新参数化不改变长度),直接用它做变分在技术上不方便。改用**能量泛函**

\[E[\gamma] = \frac{1}{2}\int_0^1 g_{ij}(\gamma)\dot{x}^i\dot{x}^j\,dt\]

它与长度泛函的关系是:能量的极值曲线自动等速参数化,且此时 \(E = \frac{1}{2}\mathcal{L}^2\),极值一致。

为什么能量泛函比长度泛函更好用? 长度泛函在根号内含有 \(\dot{x}\),使得 EL 方程非常复杂(分母有根号)。更致命的是,长度泛函具有参数化不变性——这意味着它的 EL 方程是退化的(有无穷多解,对应同一条曲线的不同参数化)。能量泛函消除了参数化自由度(它的极值曲线自动是等速参数化的),同时给出更简洁的 EL 方程。

EL 方程:对 \(L = \frac{1}{2}g_{ij}\dot{x}^i\dot{x}^j\) 写 EL 方程(对每个坐标 \(x^k\)):

\[\frac{\partial L}{\partial x^k} - \frac{d}{dt}\frac{\partial L}{\partial \dot{x}^k} = 0\]

经过 Christoffel 符号的标准计算(利用 \(\Gamma^k_{ij} = \frac{1}{2}g^{kl}(\partial_i g_{jl} + \partial_j g_{il} - \partial_l g_{ij})\)),得到

\[\boxed{\ddot{x}^k + \Gamma^k_{ij}\dot{x}^i\dot{x}^j = 0}\]

这就是**测地线方程**。在平坦空间(\(\Gamma = 0\))中退化为 \(\ddot{x}^k = 0\)(直线)。

机器人应用:SO(3) 上的测地线给出旋转的最短路径——这正是 SLERP(球面线性插值)的理论基础。SE(3) 上的测地线用于机械臂末端位姿的平滑插值。在构型空间 \(\mathcal{Q}\) 上,带度量 \(g_{ij} = M_{ij}(q)\)(惯性矩阵)的测地线是自由运动(无外力)的轨迹——这揭示了"惯性矩阵是构型空间度量"的深刻几何意义。

1.3.4 最短路径问题的统一视角

在深入最小旋转面之前,让我们从统一的视角审视已经求解的三个问题。

问题 Lagrangian \(L\) \(L\)\(x\) 简化方法 解的类型
最速降线 \(\sqrt{\frac{1+y'^2}{2gy}}\) Beltrami → 参数化 摆线
悬链线 \(y\sqrt{1+y'^2}\) Beltrami → 分离变量 双曲余弦
测地线 \(\frac{1}{2}g_{ij}\dot{x}^i\dot{x}^j\) 取决于度量 EL → Christoffel 测地线方程

所有三个问题的共同点是:通过 EL 方程或其简化形式(Beltrami),将无穷维优化问题转化为 ODE,然后通过巧妙的参数化或变量替换求解。这就是变分法的基本工作流程:

\[\text{泛函极值问题} \xrightarrow{\delta J = 0} \text{EL 方程} \xrightarrow{\text{对称性/第一积分}} \text{降阶 ODE} \xrightarrow{\text{积分}} \text{解析解}\]

需要注意的是:并非所有变分问题都有解析解。上述经典问题之所以可解,是因为它们具有特殊的对称性(\(L\) 不含 \(x\))或特殊的函数形式。对一般的 Lagrangian,EL 方程只能数值求解——这正是 1.8 节直接法存在的意义。

1.3.5 最小旋转面(Catenoid)

问题陈述:连接两个同轴圆环(半径 \(r_1\)\(r_2\),间距 \(d\))的旋转面中,哪一个面积最小?

建模:设旋转面由曲线 \(y(x)\)\(x\) 轴旋转生成。面积为

\[A[y] = 2\pi\int_a^b y\sqrt{1+y'^2}\,dx\]

注意这和悬链线问题的泛函形式相同(只差常数 \(2\pi\))!因此 EL 方程相同,解也是悬链线:

\[y(x) = c\cosh\frac{x-x_0}{c}\]

但有一个关键区别:最小旋转面问题可能没有光滑解!当两个圆环距离太远时,悬链线解不存在——这时最小面积的"解"是两个圆盘(Goldschmidt 不连续解)。这是变分法中"存在性问题"的经典例子:EL 方程有解不等于变分问题有解。

Goldschmidt 不连续解的存在条件:当 \(d/(r_1+r_2)\) 超过临界值约 0.6627 时,连续的悬链线解不存在,只有不连续的两圆盘解。

为什么这个例子在理论上如此重要? 它告诉我们变分问题的三个微妙之处:

  1. 存在性不保证:EL 方程的解(光滑极值曲线)可能不存在
  2. 正则性不保证:真正的最优解可能不是光滑的(两个圆盘是"不连续解")
  3. 全局 vs 局部:即使 EL 方程有解且满足二阶条件(局部极小),全局最小也可能是另一种完全不同的对象

这些问题催生了现代变分法中的**存在性理论**(Tonelli 直接法)、正则性理论**和**松弛理论——远超经典 EL 框架的范畴。

本质洞察:变分问题的"解"可能根本不在我们最初搜索的函数空间中。最小旋转面问题教会我们谦卑——EL 方程只是故事的一部分,完整的理论还需要存在性、正则性和全局性的分析。

1.3.6 等时降线的惊人性质(Tautochrone)

最速降线还有一个令人惊叹的推论。考虑倒置的摆线(开口向上的拱),一个质点从拱上任意高度无初速释放,滑到最低点所需的时间**与释放高度无关**——这就是等时降线性质(tautochrone)。

证明思路:利用能量守恒和摆线的参数方程,可以算出从高度 \(y_0 = r(1-\cos\alpha)\) 滑到最低点 \(\theta = \pi\) 的时间为

\[T = \pi\sqrt{\frac{r}{g}}\]

它不依赖于 \(\alpha\)(释放角度)!

工程意义:Huygens (1673) 正是利用这个性质设计了**摆线摆钟**——通过让摆锤沿摆线轨迹运动(而非圆弧),使振动周期不依赖于振幅,从而制造出精度更高的钟表。这是变分法结果第一次直接产生工程应用。

⚠️ 常见陷阱

陷阱 4(方法选择错误):对不显含 \(x\) 的问题仍然展开完整的二阶 EL 方程

  • 后果:计算量增大数倍,且容易计算出错
  • 正确做法:先检查 \(L\) 是否不含 \(x\),如果是,直接使用 Beltrami 恒等式降为一阶
  • 经验法则:看到 \(L(y, y')\)(无 \(x\))就用 Beltrami;看到 \(L(x, y')\)(无 \(y\))就用 \(L_{y'} = C\)

陷阱 5(几何直觉缺失):不理解能量泛函为何可以替代长度泛函

  • 根本原因:长度泛函有参数化自由度,能量泛函通过"固定参数化为等速"消除了这个自由度
  • 关键关系:\(E[\gamma] \geq \frac{1}{2}\mathcal{L}[\gamma]^2\)(Cauchy-Schwarz),等号当且仅当等速参数化
  • 类比:这类似于"归一化"——长度不关心你走多快,能量关心,因此能量的极值曲线同时确定了形状和参数化

陷阱 6(存在性问题):假设 EL 方程有解就意味着变分问题有解

  • 实际情况:EL 方程可能无解(如最小旋转面问题距离过大时)、有解但不是极值(鞍点)、或有多个解
  • 正确态度:EL 方程是必要条件,其解是"候选极值函数"。需要额外验证:(1) 解的存在性 (2) 二阶条件 (3) 全局比较

练习

  1. 验证摆线解确实满足 Beltrami 恒等式 \(y(1+y'^2) = 2r\)。在数值上,用 Python 画出 \(r=1\) 时的摆线并验证曲率。
  2. 对球面 \(S^2\)\(g = d\theta^2 + \sin^2\theta\,d\varphi^2\)),写出测地线方程,验证大圆弧是解。
  3. 工程联系题:minimum-snap 问题 \(J = \int_0^T |\ddddot{x}|^2\,dt\) 的 EL 方程是什么阶的 ODE?其解是什么类型的函数?解释为什么四旋翼轨迹用多项式参数化。
  4. 对最小旋转面问题,推导 Goldschmidt 条件:当两个等半径圆环 \(r_1 = r_2 = R\) 间距 \(d\) 满足什么条件时,悬链线解存在?

1.4 带约束的变分问题 ⭐⭐

动机

在实际问题中,极值函数往往需要满足额外的约束条件——不仅仅是端点约束,还有积分约束(如固定总长度)或逐点约束(如保持在某个区域内)。这正如有限维优化中的约束优化:我们需要 Lagrange 乘子法的无穷维版本。

如果不用乘子法会怎样

考虑经典的 Dido 问题:用固定长度 \(\ell\) 的绳索围出最大面积。如果不用 Lagrange 乘子法,你必须在"长度恰好为 \(\ell\) 的所有曲线"这个约束集中直接搜索,这是一个无穷维约束流形上的优化——技术上极其困难。乘子法的精妙之处在于:通过引入一个额外的标量参数 \(\lambda\),将约束优化转化为无约束优化,代价只是增加了一个未知数。

约束变分问题的分类

约束类型 数学形式 乘子性质 经典例题
等周约束 \(\int G(x,y,y')dx = \ell\) 常数 \(\lambda\) Dido 问题、等周不等式
Holonomic 约束 \(\phi(x, y_1, \ldots, y_n) = 0\) 函数 \(\lambda(x)\) 曲面上最短路径
非完整约束 \(\phi(x, y, y') = 0\)(含导数) 函数 \(\lambda(x)\)(特殊处理) 滚动无滑约束
不等式约束 \(g(x, y, y') \leq 0\) \(\lambda(x) \geq 0\)(互补松弛) 避障轨迹规划

1.4.1 等周问题(Isoperimetric Problem)

历史:等周问题是人类最早研究的优化问题之一。传说迦太基女王 Dido 请求努米底亚国王赠予"一张牛皮所能覆盖的土地",她将牛皮割成细条,围出了一大片土地——这就是 Dido 问题的传说起源。数学上,它问的是:固定周长下哪个封闭曲线围的面积最大?

问题的一般形式:在约束

\[K[y] = \int_a^b G(x, y, y')\,dx = \ell \quad (\text{给定常数})\]

下极小化 \(J[y] = \int_a^b L(x, y, y')\,dx\)

方法:构造**增广 Lagrangian**

\[\tilde{L} = L + \lambda G\]

其中 \(\lambda\) 是常数 Lagrange 乘子。对 \(\tilde{L}\) 写 EL 方程:

\[\frac{\partial(L + \lambda G)}{\partial y} - \frac{d}{dx}\frac{\partial(L + \lambda G)}{\partial y'} = 0\]

配合约束 \(K[y] = \ell\) 确定 \(\lambda\) 的值。

为什么这里 \(\lambda\) 是常数而不是函数? 因为等周约束 \(\int G\,dx = \ell\) 是一个**全局的**(积分)约束——它约束的是函数 \(y(x)\) 的一个全局性质(总长度),而不是逐点的性质。类比有限维:\(f(x) = 0\) 是一个约束方程,\(\lambda\) 是一个数;如果有 \(n\) 个约束,就有 \(n\) 个乘子。等周约束是"一个"约束(虽然是积分形式的),所以乘子是一个数。

完整推导:设 \(y^*\) 满足约束 \(K[y^*] = \ell\) 并使 \(J\) 极小。取容许变分 \(\eta_1, \eta_2\)(两个独立方向),考虑双参数族

\[y_{\varepsilon_1, \varepsilon_2} = y^* + \varepsilon_1 \eta_1 + \varepsilon_2 \eta_2\]

极值条件给出:在约束 \(K[y_{\varepsilon_1,\varepsilon_2}] = \ell\) 下,\(J\)\((\varepsilon_1, \varepsilon_2) = (0,0)\) 极小。用有限维 Lagrange 乘子法:

\[\frac{\partial J}{\partial \varepsilon_i}\bigg|_{0,0} + \lambda \frac{\partial K}{\partial \varepsilon_i}\bigg|_{0,0} = 0, \quad i=1,2\]

\(\delta J + \lambda \delta K = 0\),这等价于对 \(\tilde{J}[y] = J[y] + \lambda K[y]\) 的无约束变分为零。

经典例题:Dido 问题

陈述:给定固定长度 \(\ell\) 的曲线 \(y(x)\)(连接 \(x\) 轴上两点 \((0,0)\)\((a,0)\)),求使曲线与 \(x\) 轴围成的面积最大的形状。

  • 面积:\(J[y] = \int_0^a y\,dx\)(最大化 = 最小化 \(-J\)
  • 长度约束:\(K[y] = \int_0^a \sqrt{1 + y'^2}\,dx = \ell\)

增广 Lagrangian:\(\tilde{L} = -y + \lambda\sqrt{1+y'^2}\)

EL 方程:\(-1 - \frac{d}{dx}\frac{\lambda y'}{\sqrt{1+y'^2}} = 0\)

积分:\(\frac{\lambda y'}{\sqrt{1+y'^2}} = -(x - c_1)\)

解出 \(y'\)\(y' = \frac{-(x-c_1)}{\sqrt{\lambda^2 - (x-c_1)^2}}\)

再积分得 \(y(x) = \sqrt{\lambda^2 - (x - c_1)^2} + c_2\)

这是**圆弧方程** \((x-c_1)^2 + (y-c_2)^2 = \lambda^2\),半径为 \(|\lambda|\)

确定常数:由 \(y(0)=0\), \(y(a)=0\)\(c_1 = a/2\)(圆心在中点上方)。由长度约束确定 \(\lambda\)

深层意义:等周不等式 \(4\pi A \leq L^2\)(面积 \(A\) 与周长 \(L\) 的关系)等号成立当且仅当曲线是圆——这是变分法给出的经典几何结果。

乘子 \(\lambda\) 的物理解释\(\lambda = dJ^*/d\ell\)——约束松弛 \(\ell\) 单位增量时,最优目标值的变化率。对 Dido 问题:增加一小段绳索长度 \(\delta\ell\),面积增加约 \(\lambda\delta\ell\)\(\lambda\) 就是"绳索的边际面积产出"。这个解释在经济学(影子价格)和最优控制(共态变量的边际值解释)中有完全对应的概念。

1.4.2 Holonomic 约束

问题:当未知函数满足逐点约束 \(\phi(x, y_1(x), \ldots, y_n(x)) = 0\) 时,EL 方程如何修改?

此时 Lagrange 乘子变为函数 \(\lambda(x)\)(而非常数),修正的 EL 方程为:

\[\frac{\partial L}{\partial y_i} - \frac{d}{dx}\frac{\partial L}{\partial y_i'} + \lambda(x)\frac{\partial \phi}{\partial y_i} = 0, \quad i = 1, \ldots, n\]

\(n\) 个方程加上约束方程 \(\phi = 0\),共 \(n+1\) 个方程确定 \(n\) 个函数 \(y_i(x)\) 和乘子 \(\lambda(x)\)

为什么这里 \(\lambda\) 是函数? 因为约束 \(\phi(x,y) = 0\) 是逐点约束——在每个 \(x\) 处都必须满足。类比有限维:如果约束数等于问题维度,约束力在每一点都可以不同。\(\lambda(x)\) 就是"在 \(x\) 处维持约束所需的力"。

本质洞察:等周约束中 \(\lambda\) 是常数,因为约束 \(\int G\,dx = \ell\) 是"全局的"(积分约束);holonomic 约束中 \(\lambda(x)\) 是函数,因为约束 \(\phi(x,y) = 0\) 是"逐点的"。这个区别正是 PMP 中共态变量的前身——PMP 的共态 \(\lambda(t)\) 是满足特定 ODE 的函数,对应动力学约束 \(\dot{x} = f(x,u)\) 这个"逐时刻"的约束。

例题:曲面上的最短路径

求曲面 \(z = f(x,y)\) 上两点间的最短曲线。约束 \(\phi: z(t) - f(x(t), y(t)) = 0\)

路径长度:\(\int \sqrt{\dot{x}^2 + \dot{y}^2 + \dot{z}^2}\,dt\)

代入约束 \(\dot{z} = f_x\dot{x} + f_y\dot{y}\),得

\[\int \sqrt{(1+f_x^2)\dot{x}^2 + 2f_xf_y\dot{x}\dot{y} + (1+f_y^2)\dot{y}^2}\,dt\]

这正是度量 \(g_{ij} = \delta_{ij} + \partial_i f \partial_j f\) 下的能量泛函——约束的效果是诱导曲面上的内蕴度量,然后问题化为测地线问题。

1.4.3 多重等周约束

实际问题中常常有多个积分约束同时存在。例如,设计一座桥梁拱形:

  • 极小化弯曲能 \(J[y] = \int_a^b (y'')^2\,dx\)
  • 约束 1(固定长度):\(\int_a^b \sqrt{1+y'^2}\,dx = \ell_1\)
  • 约束 2(固定面积/载荷):\(\int_a^b y\,dx = \ell_2\)

此时需要引入**多个**乘子 \(\lambda_1, \lambda_2\)

\[\tilde{L} = L + \lambda_1 G_1 + \lambda_2 G_2\]

EL 方程加上两个约束方程,确定 \(y(x)\)\(\lambda_1\)\(\lambda_2\)

推广到 \(m\) 个约束:如果有 \(m\) 个等周约束 \(K_j[y] = \ell_j\)\(j = 1, \ldots, m\)),增广 Lagrangian 为

\[\tilde{L} = L + \sum_{j=1}^m \lambda_j G_j\]

需要满足**正则条件**:\(m\) 个变分 \(\delta K_1, \ldots, \delta K_m\) 在极值函数处**线性无关**。如果线性相关(约束退化),标准乘子法可能失效,需要考虑"异常乘子"情形(\(\lambda_0 = 0\))。

与有限维的完全对应:这和有限维约束优化中"约束规格"(LICQ/MFCQ) 的概念完全对应——约束梯度线性无关是 KKT 条件成立的前提。

1.4.5 非完整约束与变分的微妙性

问题:当约束涉及速度(如 \(\dot{y} = f(x, y)\)),情况变得微妙。

关键区别: - 完整约束 \(\phi(x,y) = 0\):约束面在构型空间中,变分 \(\delta y\) 必须切于约束面 - 非完整约束 \(a_i(q)\dot{q}_i = 0\):约束在速度空间中,\(\delta q\)\(\dot{q}\) 的约束条件**不等价**

这导致了一个深刻的区别:非完整约束下,d'Alembert 原理和 Hamilton 原理给出**不同的方程**!d'Alembert 原理(虚功原理 \(\delta q\) 满足约束)给出正确的运动方程,而 Hamilton 原理(使用 Lagrange 乘子强行约束 \(\dot{q}\))给出的是所谓的 "vakonomic" 方程——一般与物理不符。

具体例子:考虑轮式机器人的无侧滑约束 \(\dot{x}\sin\theta - \dot{y}\cos\theta = 0\)。这个约束不可积(不能化为 \(\phi(x,y,\theta) = 0\) 的形式),因此是非完整的。正确的运动方程必须用 d'Alembert 原理推导:

\[M\ddot{q} + C\dot{q} + g = B\tau + A^\top \lambda\]

其中 \(A(q)\dot{q} = 0\) 是约束方程。直接用 Hamilton 原理(将约束代入 Lagrangian 后做无约束变分)会得到错误的方程。

这个陷阱在机器人学中尤其重要:轮式移动机器人、蛇形机器人、滚动接触的操作手都涉及非完整约束。使用错误的变分原理会导致仿真与实际不符——这是 1990 年代之前文献中多次出现的错误。

如何判断约束是否完整? 检查 Frobenius 条件:约束 \(a_i(q)dq_i = 0\) 是完整的当且仅当 \(a \wedge da = 0\)(微分形式的可积条件)。如果 Frobenius 条件不满足,约束是非完整的,必须使用 d'Alembert 原理。

1.4.5 横截条件:端点约束的推广

当变分问题的端点不是完全固定时,需要额外的边界条件。以下是三种重要情形:

情形 1:一个端点自由

\(y(b)\) 自由(不固定),则变分 \(\eta\) 不再需要满足 \(\eta(b) = 0\)。分部积分后:

\[\delta J = \int_a^b \left[L_y - \frac{d}{dx}L_{y'}\right]\eta\,dx + L_{y'}\big|_{x=b} \cdot \eta(b) = 0\]

由于 \(\eta(b)\) 可以任取,而积分部分(EL 方程)已经为零,我们得到:

\[\left.\frac{\partial L}{\partial y'}\right|_{x=b} = 0 \quad (\text{自然边界条件})\]

情形 2:端点在给定曲线上

若右端点约束在曲线 \(y = \gamma(x)\) 上,即 \(y(b) = \gamma(b)\),但 \(b\) 本身可变(端点可以沿曲线滑动),则横截条件为

\[\left[L + (\gamma' - y')L_{y'}\right]_{x=b} = 0\]

情形 2 的完整推导

设端点轨迹为 \((b+\varepsilon\delta b, \gamma(b+\varepsilon\delta b))\)。第一变分包含两部分:(a) 区间内部的变化 \(\eta\);(b) 端点位置的变化 \(\delta b\)

完整的第一变分为

\[\delta J = \int_a^b [L_y - (L_{y'})'] \eta\,dx + \left[L\,\delta b + L_{y'}(\delta y - y'\delta b)\right]_{x=b}\]

其中 \(\delta y = \gamma'(b)\delta b\)(端点在曲线 \(\gamma\) 上滑动)。EL 方程消去积分部分后:

\[L(b)\,\delta b + L_{y'}(b)[\gamma'(b) - y'(b)]\delta b = 0\]

\(\delta b\) 任意性:

\[L + (\gamma' - y')L_{y'}\bigg|_{x=b} = 0\]

几何意义:如果极值曲线与约束曲线 \(\gamma\) 在端点的切线方向相同(\(y' = \gamma'\)),则横截条件退化为 \(L = 0\)。如果 \(\gamma\) 是竖直线(\(\gamma' = \infty\)),横截条件退化为 \(L_{y'} = 0\)(自然边界条件)。

情形 3:有终端代价

若问题含终端代价 \(\Phi(y(b))\):极小化 \(J[y] = \Phi(y(b)) + \int_a^b L\,dx\),则自然边界条件变为

\[L_{y'}\big|_{x=b} = -\frac{d\Phi}{dy}\bigg|_{y(b)}\]

这正是 PMP 中横截条件 \(\lambda(t_f) = \nabla\Phi(x(t_f))\) 的直接前身。

通向 PMP 的桥梁

变分法端点条件 PMP 对应 共态条件
\(y(b) = B\)(固定) \(x(t_f) = x_1\)(全固定) \(\lambda(t_f)\) 无约束(自由)
\(y(b)\) 自由 \(x(t_f)\) 自由,无终端代价 \(\lambda(t_f) = 0\)
\(y(b)\) 自由,有终端代价 \(\Phi\) \(x(t_f)\) 自由,有终端代价 \(\lambda(t_f) = \nabla\Phi(x(t_f))\)
\(y(b) \in \Gamma\) \(x(t_f) \in M_T\)(终端流形) \(\lambda(t_f) \perp T_{x(t_f)}M_T\)

本质洞察:变分法的横截条件本质上是"边界处的力平衡"——极值函数在自由端点处不能有"剩余力"推动它改变。这个力平衡条件在 PMP 中演化为共态的终端条件,在 MPC 中演化为终端约束或终端代价的选择——三者的数学本质完全相同。

1.4.7 Weierstrass-Erdmann 角条件

当极值函数不是 \(C^2\) 而只是分段 \(C^2\)(存在"角点",即导数 \(y'\) 在某些点不连续)时,EL 方程在角点处不适用。取而代之的是 Weierstrass-Erdmann 角条件

定理:设极值函数 \(y^*(x)\)\(x = c\) 处有角点(\(y'^-\)\(y'^+\) 分别为左右极限),则:

\[L_{y'}(c, y(c), y'^-) = L_{y'}(c, y(c), y'^+) \quad \text{(动量连续)}\]
\[\left[L - y'L_{y'}\right]^- = \left[L - y'L_{y'}\right]^+ \quad \text{(Hamiltonian 连续)}\]

推导:角点处 \(y\) 连续但 \(y'\) 不连续。对分段 \(C^2\) 函数,分段应用分部积分:

\[\delta J = \int_a^c [\cdots]\eta\,dx + \int_c^b [\cdots]\eta\,dx + [L_{y'}]^-_c \cdot \eta(c) - [L_{y'}]^+_c \cdot \eta(c) = 0\]

\(\eta(c)\) 的任意性得 \([L_{y'}]^- = [L_{y'}]^+\)。第二个条件(Hamiltonian 连续)来自 \(\eta\) 也可以在 \(c\) 处引起"水平方向"的变分。

物理意义:在光学中,角条件对应 Snell 折射定律——光线在界面处折射时,\(n\sin\theta\)(对应 \(L_{y'}\))连续。变分法的角条件就是 Snell 定律的推广!

最优控制中的对应:角条件对应 PMP 中的**切换条件**——当最优控制从一种模式切换到另一种模式时(如 bang-bang 控制的 \(u = +1 \to u = -1\)),共态 \(\lambda(t)\) 必须连续,Hamiltonian 也必须连续。

⚠️ 常见陷阱

陷阱 7(等周问题的异常情形):忘记检查正则条件

  • 问题:并非所有等周问题都能用 \(\tilde{L} = L + \lambda G\) 求解
  • 根本原因:如果约束 \(K[y] = \ell\) 的变分 \(\delta K\) 在极值函数处恒为零(约束退化),标准乘子法失效
  • 正确做法:先验证 \(\delta K \neq 0\)(即存在某个 \(\eta\) 使 \(\int G_y\eta + G_{y'}\eta' \neq 0\)),否则需要使用"异常乘子"理论

陷阱 8(混淆常数乘子与函数乘子)

  • 错误:对等周约束用函数 \(\lambda(x)\);或对逐点约束用常数 \(\lambda\)
  • 判断规则:积分约束→常数乘子;逐点约束→函数乘子

练习

  1. 用等周方法证明:固定周长下最大面积的封闭曲线是圆。
  2. 对曲面 \(z = f(x,y)\) 上从 \(A\)\(B\) 的最短路径问题,写出带约束 \(z - f(x,y) = 0\) 的 EL 方程系统。
  3. 跨章连接:解释等周问题中的常数乘子 \(\lambda\) 与 PMP 中的共态变量 \(\lambda(t)\) 之间的概念联系。在什么意义下,PMP 共态是"连续版本的 Lagrange 乘子"?
  4. \(J = \int_0^1 y'^2\,dx\)\(y(0)=0\)\(y(1)\) 自由,求解并解释为什么结果是 \(y \equiv 0\)(自然边界条件 \(y'(1) = 0\) 加上 EL 方程 \(y'' = 0\) 的结果)。

1.5 Hamilton 原理与经典力学 ⭐⭐

动机

变分法最伟大的应用之一是经典力学的 Lagrange 与 Hamilton 形式。这不仅在物理学中有根本性意义,更是机器人动力学方程 \(M(q)\ddot{q} + C(q,\dot{q})\dot{q} + g(q) = \tau\) 的直接来源。理解力学的变分本质,才能真正掌握机器人动力学——而非只会调用 pinocchio.rnea() 而不知其理。

Hamilton 原理的精确陈述

定理 1.5.1(Hamilton 原理 / 最小作用量原理):对 \(n\) 自由度力学系统,取广义坐标 \(q = (q_1, \ldots, q_n)\),定义 Lagrangian

\[L(q, \dot{q}, t) = T(q, \dot{q}) - V(q)\]

则**真实运动轨迹** \(q^*(t)\) 使**作用量**

\[S[q] = \int_{t_1}^{t_2} L(q, \dot{q}, t)\,dt\]

取**驻值**(\(\delta S = 0\)),即满足 EL 方程。

关键辨析:Hamilton 原理说的是"驻值"(stationary),不是"最小值"(minimum)!对短时间区间内的轨迹,作用量确实取极小值;但对长时间区间(超过第一个共轭点),可能只是鞍点。正确的表述是"驻值原理"而非"最小作用量原理"——后者虽然历史上沿用,但数学上不准确。

跨领域类比:Hamilton 原理与 Fermat 原理(光程驻值原理)的关系:

维度 Fermat 原理 Hamilton 原理
驻值量 光程 \(\int n\,ds\) 作用量 \(\int L\,dt\)
EL 方程 Eikonal/光线方程 Newton/Lagrange 方程
物理意义 光沿"最短光程"路径传播 质点沿"驻值作用量"轨迹运动
量子对应 波动光学 → 几何光学 量子力学 → 经典力学

Johann Bernoulli 正是利用 Fermat 原理和 Snell 定律的类比来解决最速降线问题的——他把重力场中的质点运动类比为折射率梯度介质中的光传播。

从 Hamilton 原理到 Lagrange 方程

\(L(q, \dot{q}, t)\) 写 EL 方程(\(q_i\) 类比 \(y\)\(t\) 类比 \(x\)):

\[\frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} - \frac{\partial L}{\partial q_i} = 0, \quad i = 1, \ldots, n\]

加上非保守力 \(Q_i\)(如摩擦、外力):

\[\frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} - \frac{\partial L}{\partial q_i} = Q_i\]

对机器人,\(Q_i = \tau_i\)(关节力矩),展开后得到标准动力学方程:

\[M(q)\ddot{q} + C(q, \dot{q})\dot{q} + g(q) = \tau\]

其中 \(M_{ij} = \frac{\partial^2 T}{\partial \dot{q}_i \partial \dot{q}_j}\) 是惯性矩阵,\(C\) 包含 Coriolis 和离心项(其元素恰好是 \(M(q)\) 在构型空间上的 Christoffel 符号),\(g = \frac{\partial V}{\partial q}\) 是重力项。

反事实推理:如果不用 Lagrange 方法而用 Newton 方法推导多体动力学会怎样?对于 \(n\) 个刚体的链式系统,Newton 方法需要为每个连接处引入约束力、对每个刚体分别列方程、然后消去约束力——方程数和工作量随 \(n\) 爆炸增长。Lagrange 方法直接在广义坐标中工作,约束力自动消失,方程数恰好等于自由度数 \(n\)。这就是为什么 Pinocchio、RBDL 等机器人动力学库本质上都基于 Lagrange/Hamilton 框架。

从 EL 方程到标准形式的完整展开

设动能 \(T = \frac{1}{2}\dot{q}^\top M(q)\dot{q}\),则

\[\frac{\partial L}{\partial \dot{q}_i} = \frac{\partial T}{\partial \dot{q}_i} = \sum_j M_{ij}(q)\dot{q}_j = [M(q)\dot{q}]_i\]
\[\frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} = \sum_j M_{ij}\ddot{q}_j + \sum_j \dot{M}_{ij}\dot{q}_j = [M\ddot{q}]_i + [\dot{M}\dot{q}]_i\]
\[\frac{\partial L}{\partial q_i} = \frac{\partial T}{\partial q_i} - \frac{\partial V}{\partial q_i} = \frac{1}{2}\sum_{j,k}\frac{\partial M_{jk}}{\partial q_i}\dot{q}_j\dot{q}_k - g_i\]

代入 EL 方程并整理,Coriolis 矩阵 \(C\) 的元素为

\[C_{ij} = \sum_k \frac{1}{2}\left(\frac{\partial M_{ij}}{\partial q_k} + \frac{\partial M_{ik}}{\partial q_j} - \frac{\partial M_{jk}}{\partial q_i}\right)\dot{q}_k\]

这恰好是 Christoffel 符号 \(\Gamma^i_{jk}\) 乘以速度 \(\dot{q}_k\) 的求和——揭示了机器人动力学方程的深层几何含义:自由运动(\(\tau = 0\), \(V = 0\))的轨迹是构型空间上以惯性矩阵为度量的测地线

Legendre 变换与 Hamilton 力学

动机:在 Lagrange 框架中,状态变量是 \((q, \dot{q})\)——位置和速度。但在最优控制中,更自然的变量是 \((q, p)\)——位置和动量。从 \((q, \dot{q})\)\((q, p)\) 的转换是通过 Legendre 变换 实现的。

Legendre 变换的几何意义:给定凸函数 \(f(v)\),它的 Legendre 变换 \(f^*(p) = \sup_v [pv - f(v)]\) 将函数"从用自变量描述"转换为"用斜率描述"。几何上:\(f(v)\) 的图形完全由其切线族确定;Legendre 变换就是用切线的斜率-截距 \((p, f^*)\) 代替点-值 \((v, f)\) 来描述同一个凸集。

定义:共轭动量 \(p_i = \frac{\partial L}{\partial \dot{q}_i}\)

Hamiltonian

\[H(q, p, t) = \sum_i p_i \dot{q}_i - L(q, \dot{q}, t)\]

其中右边的 \(\dot{q}\) 需要通过 \(p = \frac{\partial L}{\partial \dot{q}}\) 反解为 \(\dot{q} = \dot{q}(q, p)\) 后代入。这个反解的可行性要求 Legendre 条件 \(\det\left(\frac{\partial^2 L}{\partial \dot{q}_i \partial \dot{q}_j}\right) \neq 0\)(Hessian 非退化)。

EL 方程在 \((q, p)\) 变量下变为 Hamilton 正则方程

\[\dot{q}_i = \frac{\partial H}{\partial p_i}, \quad \dot{p}_i = -\frac{\partial H}{\partial q_i}\]

这是一个 \(2n\) 维的一阶 ODE 系统(与 \(n\) 维二阶系统等价)。其优美的对称结构是辛几何的核心,也是 PMP 共态方程的直接原型。

对机器人系统\(L = \frac{1}{2}\dot{q}^\top M(q)\dot{q} - V(q)\),共轭动量 \(p = M(q)\dot{q}\)(广义动量),Hamiltonian

\[H = \frac{1}{2}p^\top M(q)^{-1}p + V(q) = T + V\]

即系统的总能量。Hamilton 方程的物理意义:\(\dot{q} = M^{-1}p\)(动量→速度)和 \(\dot{p} = \tau - C^\top\dot{q} - g\)(力矩→动量变化率)。

相空间与辛结构

Hamilton 方程可以写成紧凑形式。定义相空间坐标 \(z = (q, p) \in \mathbb{R}^{2n}\) 和辛矩阵

\[J = \begin{pmatrix} 0 & I_n \\ -I_n & 0 \end{pmatrix}\]

则 Hamilton 方程为

\[\dot{z} = J \nabla_z H\]

辛结构的核心性质:Hamilton 流保持辛形式 \(\omega = \sum_i dp_i \wedge dq_i\)。几何上,这意味着相空间中的"面积"(更精确地说是辛体积)在时间演化下守恒(Liouville 定理)。

为什么辛结构对机器人学重要?

  1. 数值积分:普通的 Runge-Kutta 方法不保持辛结构,长时间模拟会导致能量漂移。辛积分器(如隐式中点法、Stormer-Verlet 法)保持辛结构,能量误差有界不累积——这对长时间仿真至关重要。
  2. 最优控制:PMP 的共态方程 \(\dot{\lambda} = -\frac{\partial H}{\partial x}\) 恰好是 Hamilton 方程的一半——最优控制问题的数学结构天然是辛的。
  3. 变分积分器:Marsden 和 West 2001 年提出的离散 Lagrange 力学,通过离散化 Hamilton 原理(而非离散化 ODE)得到保辛的积分格式。

辛积分器 vs 普通积分器的对比

考虑简谐振子 \(H = \frac{1}{2}(p^2 + q^2)\),精确解是圆轨道。

积分器 10 周期后相轨迹 100 周期后 能量误差
显式 Euler 螺旋向外(发散) 严重发散 指数增长
RK4 接近圆,轻微漂移 可见漂移 线性增长
辛 Euler(隐式) 精确圆(轻微变形) 仍然圆 有界振荡
Stormer-Verlet 精确圆 精确圆 有界振荡

对 1000 自由度的机器人仿真 10 秒钟(需要 10000 步),能量漂移的差异可以是 \(10^{-3}\)(辛积分器)vs \(10^{+2}\)(非辛积分器)——前者可用,后者完全失效。

Poisson 括号与对称性:辛结构还提供了一个优美的代数工具——Poisson 括号:

\[\{F, G\} = \sum_i \left(\frac{\partial F}{\partial q_i}\frac{\partial G}{\partial p_i} - \frac{\partial F}{\partial p_i}\frac{\partial G}{\partial q_i}\right)\]

Hamilton 方程可以写为 \(\dot{F} = \{F, H\}\)。两个守恒量的 Poisson 括号仍是守恒量(Poisson 定理)——这给出了系统地寻找新守恒量的方法。

从 Lagrange 到 Newton 的等价性:详细展示

为了让读者真正相信 Lagrange 方程与 Newton 方程是等价的,我们对最简单的情况做完整推导。

单质点在势场中的运动\(L = \frac{1}{2}m(\dot{x}^2 + \dot{y}^2 + \dot{z}^2) - V(x,y,z)\)

\(x\) 坐标写 EL 方程:

\[\frac{d}{dt}\frac{\partial L}{\partial \dot{x}} - \frac{\partial L}{\partial x} = 0\]
\[\frac{d}{dt}(m\dot{x}) - \left(-\frac{\partial V}{\partial x}\right) = 0\]
\[m\ddot{x} = -\frac{\partial V}{\partial x} = F_x\]

这正是 Newton 第二定律!\(y\)\(z\) 坐标类似。

Lagrange 方法的真正优势不在于简单系统,而在于有约束的复杂系统。考虑一个在球面 \(x^2+y^2+z^2 = R^2\) 上运动的质点。Newton 方法需要引入法向约束力 \(N\),列三个方程加一个约束方程。Lagrange 方法直接用球坐标 \((\theta, \varphi)\) 作为广义坐标:

\[L = \frac{1}{2}mR^2(\dot{\theta}^2 + \sin^2\theta\,\dot{\varphi}^2) - V(\theta, \varphi)\]

约束自动满足,方程数等于自由度数 \(2\),无需引入约束力。

反事实推理(进一步强化):考虑一个 7 自由度机械臂。Newton 方法需要对 7 个连杆分别列力与力矩方程(每个 6 个分量),加上 6 个铰链约束(每个 5 个约束方程),总共约 42 个方程和 30 个约束力——需要消去约束力后才得到 7 个运动方程。Lagrange 方法直接在 7 个关节角中工作,写出 7 个 EL 方程——不涉及任何约束力。计算量的差异是数量级的。

⚠️ 常见陷阱

陷阱 9(概念误区):把"最小作用量原理"当作物理公理

  • 实际情况:Hamilton 原理等价于 Newton 定律——两者可互推。它不是更"基本"的,只是更"方便"的
  • 但在量子力学中(Feynman 路径积分),作用量确实获得了更深层的意义——经典轨迹是所有路径中相位干涉相长的那条

陷阱 10(数学错误):Legendre 变换不存在时仍然写 Hamiltonian

  • 根本原因:当 \(\det(L_{\dot{q}\dot{q}}) = 0\) 时,\(p = L_{\dot{q}}\) 不可逆,Hamiltonian 不存在
  • 工程场景:奇异 Lagrangian(如规范场论、某些欠驱动系统)需要 Dirac 约束分析
  • 机器人学中:标准 \(L = \frac{1}{2}\dot{q}^\top M(q)\dot{q} - V(q)\) 的 Hessian 就是惯性矩阵 \(M(q)\),正定 → Legendre 变换总存在

陷阱 11(思维陷阱):认为 Lagrange 和 Hamilton 形式"只是同一个东西的改写"

  • 虽然数学上等价,但两者的**计算优势场景**完全不同
  • Lagrange 形式:适合推导具体系统的运动方程(直接写 \(T-V\)
  • Hamilton 形式:适合分析性质(守恒量、稳定性)和数值方法(辛积分器)
  • 对最优控制:Hamilton 形式是 PMP 的自然语言

练习

  1. 对简谐振子 \(L = \frac{1}{2}m\dot{x}^2 - \frac{1}{2}kx^2\),写出 Hamilton 方程并验证其与 EL 方程等价。画出相空间轨迹。
  2. 对双摆系统,写出 Lagrangian \(L(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2)\) 和两个耦合 EL 方程。
  3. 解释为什么 PMP 中的 Hamiltonian \(H = \lambda^\top f - L\) 与经典力学的 Hamiltonian \(H = p\dot{q} - L\) 形式相同但意义不同。具体地,PMP 中 \(\lambda\) 对应什么?\(f\) 对应什么?
  4. 验证 Hamilton 方程 \(\dot{z} = J\nabla H\) 保持辛形式:\(\frac{d}{dt}(dp_i \wedge dq_i) = 0\)

1.6 第一积分与守恒量 ⭐⭐

动机

求解二阶 ODE(EL 方程)通常比求解一阶 ODE 困难得多。如果能利用问题的特殊结构将二阶方程"降阶"为一阶方程,就能大大简化求解。变分法中的"第一积分"(first integral)正是这种降阶工具——它利用 Lagrangian 的对称性自动产生守恒量。

1.6.1 Beltrami 恒等式

定理 1.6.1:若 \(L = L(y, y')\) 不显含 \(x\)(即 \(\partial L / \partial x = 0\)),则沿 EL 方程的解:

\[\boxed{L - y'\frac{\partial L}{\partial y'} = C \quad (\text{常数})}\]

证明:计算

\[\frac{d}{dx}\left(L - y'L_{y'}\right) = L_y y' + L_{y'}y'' - y''L_{y'} - y'\frac{d}{dx}L_{y'}\]
\[= y'\left(L_y - \frac{d}{dx}L_{y'}\right) + \underbrace{L_{y'}y'' - y''L_{y'}}_{=0} = y' \cdot 0 = 0\]

最后一步用了 EL 方程 \(L_y - \frac{d}{dx}L_{y'} = 0\)。因此 \(L - y'L_{y'}\) 是常数。\(\square\)

物理意义:对力学系统 \(L = T(\dot{q}) - V(q)\)\(L\) 不显含 \(t\)),Beltrami 恒等式给出

\[L - \dot{q}L_{\dot{q}} = (T-V) - \dot{q}\cdot(2T/\dot{q}) = -T-V = -(T + V) = -H\]

Hamiltonian 守恒(能量守恒)。

Beltrami 恒等式的实际价值:它将二阶 ODE 降为一阶 ODE——在很多经典问题中(最速降线、悬链线、最小旋转面),这使得问题可以通过分离变量解析求解。没有 Beltrami 恒等式,这些问题仍然可解但计算量大增。

另一种有用的第一积分:若 \(L = L(x, y')\)(不含 \(y\)),则 EL 方程给出 \(\frac{d}{dx}L_{y'} = L_y = 0\),因此

\[L_{y'} = \text{const}\]

这对应"无 \(y\)→动量守恒"。

1.6.2 动量守恒

定理 1.6.2:若 \(L\) 不含某个广义坐标 \(q_k\)(即 \(\partial L / \partial q_k = 0\),称 \(q_k\) 为**循环坐标**),则共轭动量

\[p_k = \frac{\partial L}{\partial \dot{q}_k} = \text{const}\]

证明:EL 方程直接给出 \(\frac{d}{dt}p_k = \frac{d}{dt}L_{\dot{q}_k} = L_{q_k} = 0\)\(\square\)

工程应用:中心力场中角坐标 \(\varphi\) 是循环坐标 → 角动量 \(p_\varphi = mr^2\dot{\varphi}\) 守恒 → 开普勒问题降为径向一维问题。

系统化的降阶策略

  1. 检查 \(L\) 是否不显含 \(t\) → 能量守恒(Beltrami)
  2. 检查是否有循环坐标 \(q_k\) → 对应动量守恒
  3. 每个守恒量减少一个自由度
  4. \(n\) 自由度系统如果有 \(n\) 个独立守恒量(完全可积),可以通过求积法完全求解

1.6.3 Noether 定理

Noether 定理是守恒量理论的王冠:它揭示了**连续对称性**与**守恒量**之间的一一对应关系。

历史:Emmy Noether 1918 年在 Gottingen 证明了这个定理(受 Hilbert 和 Klein 关于广义相对论中能量守恒问题的讨论启发)。这个定理被 Feza Gursey 称为"20 世纪物理学中最重要的数学定理之一"。Noether 本人在当时因为女性身份无法获得正式教职——她在 Gottingen 的课程以 Hilbert 的名义开设。

定理 1.6.3(Noether 定理,1918):设 Lagrangian \(L(q, \dot{q}, t)\) 在连续变换

\[q_i \mapsto q_i + \varepsilon\xi_i(q, t), \quad t \mapsto t + \varepsilon\tau(q, t)\]

下不变(至多差一个全导数),则沿 EL 方程的解:

\[I = \sum_i \frac{\partial L}{\partial \dot{q}_i}\xi_i + \left(L - \sum_i \frac{\partial L}{\partial \dot{q}_i}\dot{q}_i\right)\tau = \text{const}\]

证明概要\(L\) 在变换下不变意味着

\[\frac{d}{d\varepsilon}\bigg|_{\varepsilon=0} L(q + \varepsilon\xi, \dot{q} + \varepsilon\dot{\xi}, t + \varepsilon\tau) = \frac{dB}{dt}\]

(某个 \(B\) 的全导数)。展开左边,利用 EL 方程 \(L_{q_i} = \frac{d}{dt}L_{\dot{q}_i}\) 消去项后,可以将全式写为某个量的 \(\frac{d}{dt}\),该量即为守恒流 \(I\)\(\square\)

对称性 \(\xi, \tau\) 守恒量 物理含义
时间平移 \(\xi=0, \tau=1\) 能量 \(H = T+V\) 系统不随时间显式变化
空间平移 \(\xi_i=\delta_{ik}, \tau=0\) 线动量 \(p_k\) 空间均匀性
旋转不变 \(\xi = \omega \times q, \tau=0\) 角动量 \(L = q \times p\) 空间各向同性
伽利略变换 \(\xi_i = t\delta_{ik}, \tau=0\) 质心运动 \(\sum m_i q_i - t P\) 惯性参考系等价

跨领域类比:Noether 定理在物理学中的地位,类似于"不变量理论"在计算机视觉中的地位——对称性减少了问题的有效自由度。正如旋转不变特征(如 SIFT)通过消除旋转自由度简化了图像匹配,Noether 守恒量通过消除对称方向简化了运动方程的求解。

反事实推理:如果没有 Noether 定理会怎样?我们仍然可以逐个验证特定守恒量(如角动量)的存在,但这需要对每个问题分别计算和验证。Noether 定理提供了一个系统化的框架:看到对称性就自动得到守恒量——这把"发现守恒量"从一种艺术变成了一种算法。

1.6.4 Noether 定理的详细例题

例题:从旋转不变性导出角动量守恒

考虑二维中心力场中的质点:

\[L = \frac{1}{2}m(\dot{x}^2 + \dot{y}^2) - V(\sqrt{x^2+y^2})\]

Step 1:识别对称性。\(L\) 在旋转 \(R_\varepsilon\) 下不变:

\[x \to x\cos\varepsilon - y\sin\varepsilon, \quad y \to x\sin\varepsilon + y\cos\varepsilon\]

Step 2:计算生成元 \(\xi\)。对 \(\varepsilon\) 求导并令 \(\varepsilon = 0\)

\[\xi_x = \frac{d}{d\varepsilon}\bigg|_0 (x\cos\varepsilon - y\sin\varepsilon) = -y\]
\[\xi_y = \frac{d}{d\varepsilon}\bigg|_0 (x\sin\varepsilon + y\cos\varepsilon) = x\]

Step 3:代入 Noether 公式(\(\tau = 0\),纯空间变换):

\[I = p_x \xi_x + p_y \xi_y = m\dot{x}(-y) + m\dot{y}(x) = m(x\dot{y} - y\dot{x})\]

这正是**角动量** \(L_z = m(x\dot{y} - y\dot{x}) = mr^2\dot{\varphi}\)

Step 4:验证。直接计算 \(\frac{dI}{dt}\)

\[\frac{dI}{dt} = m(\dot{x}\dot{y} + x\ddot{y} - \dot{y}\dot{x} - y\ddot{x}) = m(x\ddot{y} - y\ddot{x})\]

由 Newton 方程 \(m\ddot{x} = -V'(r) \frac{x}{r}\), \(m\ddot{y} = -V'(r)\frac{y}{r}\)

\[\frac{dI}{dt} = x\left(-V'\frac{y}{r}\right) - y\left(-V'\frac{x}{r}\right) = 0\]

守恒性验证通过。

这个例题展示了 Noether 定理的工作流程:(1) 识别对称性 → (2) 计算生成元 → (3) 代入公式 → (4) 验证。对于更复杂的系统(如多体系统、场论),流程完全相同。

⚠️ 常见陷阱

陷阱 12:将 Noether 定理用于离散对称性

  • 错误:反射对称 \(q \to -q\) 也应该给出守恒量
  • 实际:Noether 定理只适用于**连续单参数**变换群。离散对称性给出选择定则而非守恒流
  • 正确:离散对称性用 Ward 恒等式(量子场论)或选择规则处理

陷阱 13:忽略"至多差一个全导数"的条件

  • 有些变换下 \(L\) 不严格不变,而是 \(L \to L + \frac{d}{dt}B\)(差一个全导数项)
  • 这时守恒量要修正为 \(I - B\)
  • 例:伽利略变换下 \(L = \frac{1}{2}m\dot{x}^2\) 变化一个全导数项

练习

  1. 对中心力场 \(V(r)\),用 Noether 定理导出角动量守恒。指出对应的对称变换。
  2. 解释为什么 Beltrami 恒等式是 Noether 定理在时间平移对称性下的特例。
  3. 在 Lagrangian \(L = \frac{1}{2}m(\dot{x}^2 + \dot{y}^2) - V(x^2 + y^2)\) 中,利用旋转对称性导出角动量守恒表达式。
  4. 一个系统同时具有时间平移对称和空间旋转对称,最多有多少个独立守恒量?这些守恒量的存在对相空间轨迹的维度有什么限制?

1.7 二阶变分与充分条件 ⭐⭐⭐

动机

EL 方程只是极值的**必要条件**——类似于有限维中 \(\nabla f = 0\) 的驻点可能是极小、极大或鞍点。要确认 EL 方程的解确实是极小值(而非极大值或鞍点),我们需要检验**二阶条件**。这在工程中至关重要:如果不检验充分条件,你可能正在执行一条"最大代价"轨迹而不自知。

历史:Euler 和 Lagrange 只给出了必要条件。Legendre (1786) 首先研究二阶条件但其论证有缺陷。Jacobi (1837) 引入共轭点理论修补了 Legendre 的不足。Weierstrass (1879) 给出了处理强极小的 E-函数方法——完成了变分法充分条件理论的大厦。整个过程经历了近一个世纪,说明充分条件远比必要条件困难。

1.7.1 第二变分的推导

定义第二变分为

\[\delta^2 J[y; h] = \frac{d^2}{d\varepsilon^2}\bigg|_{\varepsilon=0} J[y + \varepsilon h]\]

对标准泛函 \(J[y] = \int_a^b L(x,y,y')dx\),我们来详细计算。设 \(\Phi(\varepsilon) = J[y + \varepsilon h]\),则

\[\Phi'(\varepsilon) = \int_a^b [L_y h + L_{y'} h'] dx\]
\[\Phi''(\varepsilon) = \int_a^b [L_{yy}h^2 + 2L_{yy'}hh' + L_{y'y'}(h')^2] dx\]

\(\varepsilon = 0\) 处取值(所有偏导数在极值函数 \(y^*\) 处计算):

\[\delta^2 J = \int_a^b \left[P(x)h^2 + 2Q(x)hh' + R(x)(h')^2\right]dx\]

其中 \(P = L_{yy}\)\(Q = L_{yy'}\)\(R = L_{y'y'}\)(均沿极值曲线取值)。

极小的**二阶必要条件**是 \(\delta^2 J \geq 0\) 对所有满足 \(h(a) = h(b) = 0\) 的容许变分 \(h\) 成立。

类比有限维:这等价于有限维中 Hessian 矩阵半正定 \(\nabla^2 f \succeq 0\) 的条件。但无穷维中"半正定"的验证远比有限维复杂——无穷维的"Hessian"是一个算子,其正定性不能通过有限步计算确认。

1.7.2 Legendre 必要条件(强/弱形式)

定理 1.7.1(Legendre 必要条件):若 \(y^*\) 是弱极小值,则

\[\boxed{L_{y'y'}(x, y^*(x), y^{*\prime}(x)) \geq 0, \quad \forall x \in [a,b]}\]

完整证明

取高频测试函数。设 \(x_0 \in (a,b)\),取光滑 bump 函数 \(\eta(x)\) 支撑在 \((x_0-\delta, x_0+\delta)\) 内,\(\eta \geq 0\),且 \(\eta(x_0) = 1\)。定义

\[h_n(x) = \frac{1}{n}\eta(x)\sin(n(x-x_0))\]

\(h_n \to 0\) 一致收敛(\(C^0\)\(C^1\) 范数趋于零),而

\[(h_n')^2 = \eta^2\cos^2(n(x-x_0)) + O(1/n)\]

代入第二变分:

\[\delta^2 J[y^*; h_n] = \int_a^b R(x) \eta^2(x)\cos^2(n(x-x_0))\,dx + O(1/n)\]

由 Riemann-Lebesgue 引理,\(\cos^2\) 的平均值趋于 \(1/2\)

\[\lim_{n\to\infty} \delta^2 J[y^*; h_n] = \frac{1}{2}\int_{x_0-\delta}^{x_0+\delta} R(x)\eta^2(x)\,dx\]

\(R(x_0) < 0\),则对足够小的 \(\delta\),上式为负 → 矛盾。故 \(R(x) = L_{y'y'} \geq 0\)\(\square\)

Legendre 的强形式:若 \(L_{y'y'} > 0\)(严格正),称满足**强化 Legendre 条件**。强化条件保证 \(L\) 关于 \(y'\) 是严格凸的——这是后续 Jacobi 充分条件的前提。

物理意义\(L_{y'y'} \geq 0\) 意味着 Lagrangian 关于 \(y'\) 是凸的——这保证了"速度方向"上的稳定性。类比:有限维中 Hessian \(\nabla^2 f \succeq 0\) 保证极值点附近是"碗形"的。凸性意味着"斜率的增加需要更多的代价"——这正是极小的基本要求。

1.7.3 Jacobi 方程与共轭点

Legendre 条件保证了"局部"(高频方向)的非负性,但"全局"(低频方向)的非负性还需要 Jacobi 条件。

为什么 Legendre 条件不够? Legendre 条件是逐点的:它保证 \(R(x) \geq 0\) 在每个 \(x\) 处。但 \(\delta^2 J\) 是一个积分——即使 \(R \geq 0\),如果 \(P\) 项(含 \(h^2\))足够负且 \(h\) 变化缓慢(\(h'\) 小),整个积分仍可能为负。Jacobi 条件正是处理这种"低频"方向的工具。

推导 Jacobi 方程:对第二变分做分部积分,将其写为

\[\delta^2 J = \int_a^b [Ph^2 + 2Qhh' + R(h')^2] dx\]

\(Rh'\) 做分部积分:\(\int R(h')^2 dx = [Rh'h]_a^b - \int (Rh')' h\,dx\)(边界项为零)

更系统的方法是对 \(\delta^2 J\) 本身做变分(对 \(h\) 的 EL 方程),得到:

定义(Jacobi 方程 / accessory equation):沿极值曲线 \(y^*\) 定义

\[\frac{d}{dx}\left(R(x) \cdot u'\right) - \left(P(x) - Q'(x)\right) \cdot u = 0\]

\[\frac{d}{dx}\left(L_{y'y'} \cdot u'\right) - \left(L_{yy} - \frac{d}{dx}L_{yy'}\right) \cdot u = 0\]

这是关于 \(u(x)\) 的二阶线性 ODE。

定义(共轭点):Jacobi 方程的解 \(u(x)\)(满足 \(u(a) = 0\), \(u'(a) \neq 0\))在 \((a,b)\) 内首次为零的点 \(\tilde{x}\) 称为 \(a\) 的**共轭点**。

定理 1.7.2(Jacobi 必要条件):若 \(y^*\) 是弱极小值,则 \((a, b)\) 内不存在共轭点(即 Jacobi 方程的解在 \((a,b)\) 内不为零)。

定理 1.7.3(Jacobi 充分条件):若 1. \(L_{y'y'} > 0\)\([a,b]\) 上(强化 Legendre) 2. \((a, b]\) 内无共轭点

\(y^*\) 是**严格弱极小值**。

共轭点的几何直觉

考虑从 \(a\) 出发的一族极值曲线(满足 EL 方程但初始斜率不同)。这些曲线初始时发散,但可能在某个点 \(\tilde{x}\) 重新汇聚——就像凸透镜的焦点。\(\tilde{x}\) 就是共轭点。

在共轭点之前,从 \(a\) 出发的极值曲线是局部最优的——没有"邻近的竞争者"能做得更好。过了共轭点之后,存在另一条极值曲线(从 \(a\) 出发,在 \(\tilde{x}\) 处与原曲线重合),使得原曲线不再最优。

经典例子:球面上的测地线。从北极出发的所有经线(大圆弧)在南极(对径点)重新汇聚→南极是北极的共轭点。因此长度超过半个大圆的测地线不再是最短路径——从另一侧绕过去更短。

1.7.4 Weierstrass E-函数与强极小

为什么还需要 Weierstrass 条件? Legendre-Jacobi 条件保证的是**弱极小**——只在 \(C^1\) 范数意义下相邻的竞争者中最优。但如果允许竞争者具有大的导数变化(\(C^0\) 范数意义下相邻),弱极小可能不是强极小。Weierstrass 1879 年的贡献是给出处理强极小的工具。

定义:Weierstrass 超越函数(excess function)

\[\mathcal{E}(x, y, y', p) = L(x, y, p) - L(x, y, y') - (p - y')L_{y'}(x, y, y')\]

几何意义\(\mathcal{E}\) 衡量的是 \(L\) 作为 \(y'\) 的函数在点 \((x,y,y')\) 处的切线与函数本身的差。具体地:\(L(x,y,y') + (p-y')L_{y'}\)\(L\)\(y'\) 处的一阶 Taylor 展开(切线近似)。\(\mathcal{E}\) 就是真实值 \(L(x,y,p)\) 与切线值之差。

因此:\(\mathcal{E} \geq 0\) 等价于 \(L\) 关于第三个变量(斜率 \(p\))是凸的——曲线在切线之上。

与 Legendre 条件的关系:对 \(p\) 接近 \(y'\)(小偏差),Taylor 展开给出

\[\mathcal{E} \approx \frac{1}{2}L_{y'y'}(p - y')^2\]

因此 \(L_{y'y'} \geq 0\)(Legendre)等价于 \(\mathcal{E} \geq 0\) 对**小偏差** \(p-y'\)。Weierstrass 条件要求 \(\mathcal{E} \geq 0\) 对**所有** \(p\)——这是更强的条件。

定理 1.7.4(Weierstrass 充分条件):若沿极值曲线 \(y^*\): 1. \(\mathcal{E}(x, y^*, y^{*\prime}, p) \geq 0\) 对所有 \(p\) 成立 2. 极值曲线可嵌入一个"场"(Hamilton-Jacobi 理论的斜率场覆盖条件)

\(y^*\) 是**强极小值**。

充分条件总结

条件 数学表述 保证的极值类型 适用场景 检验难度
Legendre(弱) \(L_{y'y'} \geq 0\) 弱极小必要 初步筛选 容易
Legendre(强化) \(L_{y'y'} > 0\) 充分条件前提 排除退化 容易
Jacobi \((a,b)\) 内无共轭点 弱极小(配合强化 Legendre 充分) 确认弱极小 中等
Weierstrass \(\mathcal{E} \geq 0 \; \forall p\) 强极小 确认全局凸性 较难
场覆盖 极值曲线族无焦点 与 Weierstrass 配合 理论完备性 较难

实践中的检验策略

  1. 先算 \(L_{y'y'}\):如果它处处 \(>0\)(很多物理问题满足),则 \(L\) 关于 \(y'\) 严格凸 → Weierstrass \(\mathcal{E} \geq 0\) 自动满足
  2. 再解 Jacobi 方程检查共轭点:如果区间短(\(b-a\) 小),通常无共轭点
  3. 对工程问题,常常不需要严格验证充分条件——物理直觉加上数值验证(微扰后泛函值增大)通常足够

⚠️ 常见陷阱

陷阱 14:把 Legendre 条件当作充分条件

  • \(L_{y'y'} \geq 0\) 只是必要条件,不是充分条件!
  • 还需要 Jacobi 条件(无共轭点)才能确认弱极小
  • 例:\(J = \int_0^T (y'^2 - y^2)\,dx\)\(L_{y'y'} = 2 > 0\),但当 \(T > \pi\) 时有共轭点,解不再是极小

陷阱 15:混淆弱极小和强极小的充分条件

  • Legendre + Jacobi → 弱极小(只在 \(C^1\) 邻域中最优)
  • Weierstrass + 场覆盖 → 强极小(在 \(C^0\) 邻域中最优)
  • 对最优控制中的 bang-bang 控制,需要强极小条件

1.7.5 经典例题:Jacobi 条件的应用

例题:考虑 \(J[y] = \int_0^T [(y')^2 - y^2]\,dx\)\(y(0) = 0\)\(y(T) = 0\)

Step 1:EL 方程 \(y'' + y = 0\),通解 \(y = A\sin x + B\cos x\)

边界条件 \(y(0) = 0\)\(B = 0\)\(y(T) = 0\)\(A\sin T = 0\)

如果 \(T \neq n\pi\),则 \(A = 0\),唯一解 \(y^* \equiv 0\)。 如果 \(T = n\pi\),则 \(A\) 任意,有无穷多解。

Step 2:对 \(y^* \equiv 0\) 检验 Legendre 条件。\(L = (y')^2 - y^2\)\(L_{y'y'} = 2 > 0\)。满足。

Step 3:检验 Jacobi 条件。沿 \(y^* = 0\)\(P = L_{yy} = -2\)\(R = L_{y'y'} = 2\)\(Q = L_{yy'} = 0\)

Jacobi 方程:\((2u')' - (-2)u = 0\),即 \(u'' + u = 0\)

\(u(x) = \sin x\)(满足 \(u(0) = 0\), \(u'(0) = 1\))。

\(u(x)\)\(x = \pi\) 处首次为零 → \(\pi\) 是共轭点。

结论: - \(T < \pi\):无共轭点,\(y^* = 0\) 是严格弱极小(Legendre + Jacobi 充分) - \(T = \pi\):共轭点恰在端点,临界情况 - \(T > \pi\)\((0, T)\) 内存在共轭点 \(\pi\)\(y^* = 0\) **不是**极小值

物理验证:当 \(T > \pi\) 时,取 \(h = \sin(\pi x/T)\),计算第二变分:

\[\delta^2 J = \int_0^T [2(\pi/T)^2\cos^2(\pi x/T) - 2\sin^2(\pi x/T)]\,dx = T[(\pi/T)^2 - 1]\]

\(T > \pi\)\((\pi/T)^2 < 1\)\(\delta^2 J < 0\) → 确实不是极小。

练习

  1. 对最速降线问题验证 Legendre 条件 \(L_{y'y'} \geq 0\) 是否满足。
  2. 解释为什么球面上长度大于半周的大圆弧不是最短路径(用共轭点的语言)。
  3. \(J = \int_0^T (y'^2 - y^2)\,dx\)\(y(0)=0\), \(y(T)=0\),用 Jacobi 条件确定 \(T\) 的临界值(在什么 \(T\) 值时出现第一个共轭点)。
  4. \(J = \int_0^1 [(y')^2 + \alpha y^2]\,dx\)\(\alpha\) 为参数),讨论当 \(\alpha\) 取不同值时 Jacobi 条件的变化。临界值 \(\alpha_c\) 是什么?

1.8 直接法与数值变分 ⭐⭐⭐

动机

到目前为止,我们使用的都是**间接法**:先写出 EL 方程(必要条件),再求解 ODE。这种方法对简单问题(一维、特殊结构)很有效,但对复杂工程问题有严重局限:

  1. EL 方程可能无解析解
  2. BVP 的数值求解可能不收敛(对初值敏感)
  3. 不等式约束很难在 EL 框架中处理

直接法(direct method)提供了另一条路径:不写 EL 方程,而是直接在有限维子空间中近似最优函数。这是现代工程计算(有限元法、轨迹优化)的数学基础。

1.8.1 Ritz 法

核心思想(Walter Ritz, 1908):选择一组基函数 \(\{\varphi_k(x)\}_{k=1}^N\),将未知函数近似为

\[y_N(x) = y_0(x) + \sum_{k=1}^N c_k \varphi_k(x)\]

其中 \(y_0(x)\) 满足边界条件,\(\varphi_k(a) = \varphi_k(b) = 0\)。代入泛函后,

\[J[y_N] = J(c_1, \ldots, c_N) \equiv \Phi(c_1, \ldots, c_N)\]

变为 \(N\) 个实数的普通函数。对 \(\Phi\) 求极值:

\[\frac{\partial \Phi}{\partial c_k} = 0, \quad k = 1, \ldots, N\]

得到 \(N\) 个代数方程确定 \(c_1, \ldots, c_N\)

Ritz 法的关键选择

选择 要求 常用选项 优劣
基函数 \(\varphi_k\) 满足边界条件、线性无关 三角函数、多项式、B-样条 三角函数收敛快(Fourier 理论);多项式简单但高阶不稳定
基函数数目 \(N\) 越大越精确,但计算量 \(O(N^2)\) 工程中 \(N = 5\sim 50\) 太少不准确,太多有 Runge 现象

例题:用 Ritz 法近似求解 \(J[y] = \int_0^1 [(y')^2 + y^2 - 2xy]\,dx\)\(y(0)=y(1)=0\)

\(\varphi_1(x) = x(1-x)\)(满足边界条件的最简多项式),\(y_1 = c_1 x(1-x)\)

\[J[c_1 x(1-x)] = \int_0^1 [c_1^2(1-2x)^2 + c_1^2 x^2(1-x)^2 - 2xc_1 x(1-x)]dx\]

计算各积分并对 \(c_1\) 求导令其为零,得到 \(c_1\) 的值。

完整的 Ritz 法计算过程(以 \(N=2\) 为例):

\(\varphi_1(x) = \sin(\pi x)\)\(\varphi_2(x) = \sin(2\pi x)\)\(y_0 = 0\)(零边界条件)。

\[y_2(x) = c_1\sin(\pi x) + c_2\sin(2\pi x)\]

代入 \(J = \int_0^1 [(y')^2 + y^2 - 2xy]\,dx\)

\[J(c_1, c_2) = \int_0^1 [(c_1\pi\cos\pi x + c_2 2\pi\cos 2\pi x)^2 + (c_1\sin\pi x + c_2\sin 2\pi x)^2 - 2x(c_1\sin\pi x + c_2\sin 2\pi x)]\,dx\]

利用三角函数的正交性 \(\int_0^1 \sin(m\pi x)\sin(n\pi x)\,dx = \frac{1}{2}\delta_{mn}\),得

\[J(c_1,c_2) = \frac{1}{2}(\pi^2+1)c_1^2 + \frac{1}{2}(4\pi^2+1)c_2^2 - c_1\int_0^1 2x\sin(\pi x)\,dx - c_2\int_0^1 2x\sin(2\pi x)\,dx\]

\(c_1, c_2\) 分别求导令其为零:

\[(\pi^2+1)c_1 = 2\int_0^1 x\sin(\pi x)\,dx = \frac{2}{\pi}\]
\[(4\pi^2+1)c_2 = 2\int_0^1 x\sin(2\pi x)\,dx = -\frac{1}{\pi}\]

解得 \(c_1 = \frac{2}{\pi(\pi^2+1)} \approx 0.0585\)\(c_2 = \frac{-1}{\pi(4\pi^2+1)} \approx -0.00806\)

这个计算过程展示了 Ritz 法的核心:泛函极值问题被完全转化为多元函数的极值问题——后者我们在大一微积分中就已经掌握。

收敛性定理:如果基函数系在函数空间中是完备的(即有限线性组合可以逼近任何容许函数),则 Ritz 近似 \(y_N\)\(N \to \infty\) 收敛到真实极值函数。这是 Sobolev 空间理论保证的。

收敛速度:收敛速度取决于解的光滑性和基函数的选择。对光滑解,三角函数基(Fourier 级数)给出指数收敛;多项式基给出代数收敛。对有角点或不连续导数的解,收敛会变慢(Gibbs 现象)。

1.8.2 Galerkin 法

与 Ritz 法的区别:Ritz 法直接极小化泛函;Galerkin 法则要求 EL 方程的残差与基函数正交。

设 EL 方程为 \(\mathcal{F}[y] := L_y - \frac{d}{dx}L_{y'} = 0\)。Galerkin 法要求:

\[\int_a^b \mathcal{F}[y_N](x) \cdot \varphi_k(x)\,dx = 0, \quad k = 1, \ldots, N\]

即残差在基函数张成的子空间上投影为零("加权残差法")。

Ritz vs Galerkin 的关系

  • 对线性问题(\(L\) 关于 \(y, y'\) 是二次的),两者等价
  • 对非线性问题,两者一般不等价,但都收敛到真实解
  • Galerkin 法更通用(不要求问题有变分形式)

本质洞察:Ritz 法本质上是"先离散再优化"——在有限维子空间中优化泛函。Galerkin 法本质上是"先优化再离散"——写出无穷维方程(EL),然后在有限维子空间中求近似解。这两种路径在变分法中等价,但在最优控制中导致了"直接法 vs 间接法"的根本分野。

1.8.3 有限元法的变分基础

有限元法(FEM)可以看作 Ritz/Galerkin 法的一种特殊实现:

  1. 分片基函数:不用全局基函数,而是把区间分成小段,在每段上用低阶多项式
  2. 局部支撑:每个基函数只在少数相邻单元上非零 → 刚度矩阵稀疏
  3. 收敛的两种途径:增加单元数(h-refinement)或增加多项式阶(p-refinement)

变分法到 FEM 的逻辑链

\[\text{物理问题} \xrightarrow{\text{变分原理}} \delta J = 0 \xrightarrow{\text{Galerkin}} \text{弱形式} \xrightarrow{\text{FEM}} K\mathbf{c} = \mathbf{f}\]

每一步的含义: - 变分原理:物理律等价于泛函驻值 - Galerkin:在有限维子空间中求近似 - FEM:选择分片多项式基函数,得到稀疏线性系统

机器人学中的应用:弹性体仿真(软体机器人)、结构分析、接触力学的有限元都基于变分原理。

从变分法到 FEM 的完整逻辑链

  1. 物理问题:弹性体在外力下的变形
  2. 变分原理:平衡构型使势能 \(\Pi[u]\) 取驻值(\(\delta\Pi = 0\)
  3. 弱形式\(\int_\Omega \sigma_{ij}\delta\epsilon_{ij}\,dV = \int_\Omega f_i\delta u_i\,dV + \int_{\Gamma_N} t_i\delta u_i\,dS\)(虚功原理)
  4. Galerkin 离散\(u_h = \sum_I N_I(x) u_I\)(用形函数 \(N_I\) 近似位移场)
  5. 刚度矩阵\(K_{IJ} = \int_\Omega B_I^\top D B_J\,dV\)\(B\) 是应变-位移矩阵,\(D\) 是本构矩阵)
  6. 线性系统\(K\mathbf{u} = \mathbf{f}\)

每一步都有明确的数学根据——变分原理保证了近似的最优性(最佳逼近性质),Galerkin 正交性保证了误差最小化。

1.8.4 配点法(Collocation Method)

配点法是直接法的另一种实现,它与 Ritz/Galerkin 有微妙的区别:

方法 要求 优势
Ritz 极小化泛函 \(J[y_N]\) 直接优化目标
Galerkin 残差与基函数正交 不需要变分形式
配点法 EL 方程在选定点处精确成立 实现最简单

配点法选择 \(N\) 个**配点** \(x_1, \ldots, x_N\),要求近似解 \(y_N\) 在这些点上精确满足 EL 方程:

\[\mathcal{F}[y_N](x_k) = 0, \quad k = 1, \ldots, N\]

与轨迹优化的联系:现代轨迹优化中的**直接配点法**(direct collocation),如 DIRCOL、Hermite-Simpson 配点,正是变分法配点法在最优控制中的推广。CasADi、Drake 等工具的配点法求解器,其数学根源就在这里。

配点法的工作流程

  1. 将时间/空间离散为 \(N\) 个节点
  2. 在每个节点上将状态/控制参数化为未知数
  3. 要求动力学方程在配点上精确成立(等式约束)
  4. 代价函数用求积公式近似(如 Gauss-Lobatto)
  5. 得到一个 NLP 问题,用 IPOPT/SNOPT 求解

这正是"先离散再优化"(直接法)的具体实现——与变分法"先优化再离散"(间接法)形成对偶。

1.8.5 直接法与间接法的对比

维度 间接法(EL → ODE → 解) 直接法(离散 → NLP → 解)
思路 先导出最优性条件,再求解 先离散化,再数值优化
数学基础 EL 方程、BVP Ritz/Galerkin、NLP
计算精度 可以达到任意精度 受网格/基函数限制
收敛性 BVP 对初值敏感 NLP 有成熟全局化技术
约束处理 不等式约束困难 容易添加不等式约束
共态信息 自然提供 \(\lambda(t)\) 需要从 KKT 乘子恢复
工程应用 低维解析解、航天 机器人 MPC、轨迹优化

反事实推理:如果 Ritz/Galerkin 法不存在会怎样?我们只能用间接法——写出 EL 方程,然后求解 BVP。对于复杂的多体系统(如人形机器人),EL 方程是高度非线性的耦合 ODE 系统,BVP 的数值求解极其困难(shooting 方法对初值猜测极度敏感)。直接法的出现使得"先离散再优化"成为可能——这正是现代轨迹优化(CasADi、Crocoddyl)的基础。

⚠️ 常见陷阱

陷阱 16(Ritz 法的基函数选择错误):基函数不满足边界条件

  • 错误做法:取 \(\varphi_k(x) = x^k\)(对 \(y(0) = 0, y(1) = 0\) 问题,\(\varphi_k(1) = 1 \neq 0\)
  • 正确做法:基函数必须满足齐次边界条件 \(\varphi_k(a) = \varphi_k(b) = 0\)
  • 典型安全选择:\(\sin(k\pi x)\)\(x(1-x)P_k(x)\)\(P_k\) 是多项式)

陷阱 17(数值稳定性):高阶多项式基函数导致病态

  • 问题:当 \(N\) 较大时,\(\{x^k\}\) 基函数组成的矩阵条件数指数增长
  • 根本原因:高阶单项式在 \([0,1]\) 上几乎线性相关(Runge 现象的根源)
  • 正确做法:使用正交基(Legendre 多项式、三角函数)或分片基(FEM)

练习

  1. 用 Ritz 法(取 \(\varphi_1 = \sin(\pi x)\))近似求解 \(J = \int_0^1 [(y')^2 - y^2]\,dx\)\(y(0)=y(1)=0\)。与精确解对比。
  2. 解释为什么有限元法使用分片线性基函数比全局多项式基函数在实践中更好(提示:稀疏性、局部支撑、条件数)。
  3. 跨章综合题:考虑 minimum-energy 问题 \(\min \int_0^T u^2\,dt\),约束 \(\dot{x} = u\)\(x(0) = 0\)\(x(T) = 1\)。(a) 用间接法(EL 方程)求解析解;(b) 用 Ritz 法(\(u(t) = c_0 + c_1 t\))求近似解;(c) 比较两者。
  4. \(J = \int_0^1 (y''^2)\,dx\)\(y(0)=y(1)=0\)\(y'(0)=1\)\(y'(1)=0\),用分段三次多项式(3个段)实现 Ritz 法并求解。这就是最简单的样条有限元。

1.9 工程应用与通向最优控制的桥梁 ⭐⭐⭐

动机

变分法不是博物馆里的古董——它是现代最优控制、机器人动力学、轨迹优化的数学根基。本节建立从变分法到 Pontryagin 极大值原理(PMP)的逻辑桥梁,让读者理解:PMP 不是从天而降的独立理论,而是变分法在面对新困难(控制约束)时的自然推广。

间接法 vs 直接法

变分法给出的 EL 方程/PMP 是一种**间接法**(indirect method):先推导最优性的必要条件(微分方程),再求解该方程。与之对立的是**直接法**(direct method):先把连续优化问题离散化为有限维 NLP,再用数值优化求解。

维度 间接法(变分→ODE) 直接法(离散→NLP)
思路 先优化再离散 先离散再优化
数学基础 EL 方程、PMP、BVP KKT 条件、NLP 求解器
优势 精度高、有显式共态信息 收敛域大、易处理不等式约束
劣势 共态初值敏感、不等式约束困难 精度受网格限制
典型工具 solve_bvp、shooting CasADi + IPOPT、Drake
机器人应用 低维解析解、航天 实时 MPC、全身控制

工业实践:机器人轨迹优化领域几乎 100% 使用直接法(CasADi、acados、Crocoddyl、Drake),间接法仅用于教学和低维解析解。但理解间接法(变分法/PMP)对于正确理解直接法的数学结构是不可或缺的——直接法的 KKT 乘子就是 PMP 共态的离散近似。

机器人动力学的变分基础

机器人标准动力学方程 \(M(q)\ddot{q} + C(q,\dot{q})\dot{q} + g(q) = \tau\) 的三种等价推导路径:

  1. Newton-Euler 递推:逐杆计算力与力矩,效率高但物理直觉差
  2. Lagrange 方法:从 \(L = T - V\) 出发写 EL 方程,优雅但对高自由度系统展开困难
  3. Hamilton 方法:在相空间 \((q, p)\) 中写正则方程,适合理论分析(辛积分器、PMP)

三者在数学上完全等价,但在数值实现和直觉方面各有优劣。Pinocchio/RBDL 使用 Newton-Euler 递推(RNEA 算法,\(O(n)\) 复杂度)计算逆动力学,但其理论基础仍然是 Lagrange 力学。

Minimum-snap 轨迹优化

四旋翼轨迹优化中,thrust 正比于加速度的导数(jerk),推力变化率正比于 snap(加速度的二阶导数)。为使飞行平滑,极小化

\[J = \int_0^T \|\ddddot{x}(t)\|^2\,dt\]

EL 方程(Euler-Poisson 方程)给出八阶 ODE \(x^{(8)} = 0\),解为**七次多项式**。这就是 Mellinger & Kumar 2011 工作的数学核心——minimum-snap 轨迹是分段七次多项式。

为什么是七次多项式? \(x^{(8)} = 0\) 意味着 \(x(t)\)\(t\) 的不超过 7 次的多项式。八阶 ODE 需要八个边界条件:起点和终点各提供位置、速度、加速度和 jerk 四个条件。

实际实现中的 Ritz 法联系:实际的 minimum-snap 求解器不是直接解 BVP,而是使用 QP 形式:将轨迹参数化为分段多项式(Ritz 法的精神),代入 snap 积分后得到关于多项式系数的二次代价函数,加上连续性和边界约束后变为一个 QP 问题。这正是"直接法"的思路——先离散(参数化为多项式),再优化(解 QP)。

从 minimum-snap 到一般轨迹优化的推广

优化目标 Lagrangian \(L\) EL 阶数 解的类型 应用
Minimum-jerk \(\|\dddot{x}\|^2\) 6 阶 5 次多项式 机械臂平滑运动
Minimum-snap \(\|x^{(4)}\|^2\) 8 阶 7 次多项式 四旋翼飞行
Minimum-crackle \(\|x^{(5)}\|^2\) 10 阶 9 次多项式 极端平滑需求
Minimum-energy \(\|u\|^2\)\(\ddot{x}=u\) 4 阶 3 次多项式 能量最优

所有这些问题的数学结构完全相同——只是 Lagrangian 中导数的阶数不同。变分法的统一框架使我们可以用同一套工具处理所有这些问题。

弹性力学中的变分原理

变分法在连续介质力学中同样发挥核心作用。弹性体的平衡构型使势能泛函取驻值:

\[\Pi[u] = \frac{1}{2}\int_\Omega \sigma_{ij}\epsilon_{ij}\,dV - \int_\Omega f_i u_i\,dV - \int_{\Gamma_N} t_i u_i\,dS\]

其中 \(u\) 是位移场,\(\sigma\) 是应力,\(\epsilon\) 是应变,\(f\) 是体力,\(t\) 是面力。

\(\delta\Pi = 0\) 给出平衡方程(EL 方程的多维推广):

\[\frac{\partial\sigma_{ij}}{\partial x_j} + f_i = 0\]

有限元法正是对这个变分原理的 Galerkin 近似——这解释了为什么有限元法如此自然地适用于结构分析:它直接近似变分原理,而非近似微分方程。

从 EL 到 PMP 的推广路径

变分法(EL 方程)的核心假设是:控制变量无约束——使 \(\delta J = 0\) 可以通过 \(\partial H/\partial u = 0\) 实现。但在真实工程中:

  • 推力有限:\(0 \leq u \leq u_{\max}\)
  • 力矩有限:\(|\tau_i| \leq \tau_{\max}\)
  • 甚至控制可能是离散的:\(u \in \{0, 1\}\)

当最优控制恰好在约束边界上时(如 bang-bang 控制),\(\partial H/\partial u = 0\) 无解——极值点不在内部!

PMP 的核心突破:用"在约束集 \(U\) 上极大化 Hamiltonian"

\[H(x, \lambda, u^*) = \max_{u \in U} H(x, \lambda, u)\]

替代"Hamiltonian 对 \(u\) 的梯度为零"。这个极大化条件在 \(u^*\) 位于 \(U\) 内部时退化为 \(\partial H/\partial u = 0\)(回归 EL),在边界时给出 bang-bang 或奇异控制。

从变分法到 PMP 的概念对应

变分法 PMP 区别/推广
\(y(x)\)(待优化函数) \(u(t)\)(控制函数) PMP 允许 \(u\) 有约束
EL 方程 \(L_y - (L_{y'})' = 0\) 共态方程 \(\dot{\lambda} = -H_x\) 形式相同,本质相同
Beltrami:\(H\) = const PMP:\(H\) 沿最优轨迹恒定 完全对应
自然边界条件 横截条件 \(\lambda(t_f) = \nabla\Phi\) 完全对应
\(\frac{\partial L}{\partial y'} = 0\)(驻值) \(\max_u H\)(极大化) PMP 的核心推广
Legendre 条件 \(L_{y'y'} \geq 0\) \(H_{uu} \leq 0\)(Hamiltonian 关于 \(u\) 凹) 符号相反(因为是极大化)

为什么 PMP 用"极大化"而非"极小化"? 这是历史约定。Pontryagin 团队原始定义的 Hamiltonian 是 \(H = \lambda^\top f + L\)(注意符号),对它极大化等价于对 \(-H\) 极小化。不同教材对 \(H\) 的符号约定不同,但物理本质相同。常见的两种约定如下:

约定 Hamiltonian 定义 最优性条件 使用的教材
Pontryagin 原始 \(H = \lambda^\top f - L_0\) \(\max_u H\) Pontryagin 1962
现代控制论 \(H = \lambda^\top f + L_0\) \(\min_u H\) Kirk, Liberzon

两者通过 \(\lambda \to -\lambda\) 互相转换,物理结果完全一致。初学者最常犯的错误就是混用两种约定——在一个问题中必须始终坚持同一种。

本质洞察:变分法与 PMP 的关系,如同无约束优化与约束优化的关系。变分法(\(\nabla J = 0\))对应无约束的 \(\nabla f = 0\);PMP(\(\max_u H\))对应约束集上的 KKT 条件——当最优解在内部时,KKT 退化为无约束条件;当最优解在边界时,KKT 给出额外的互补松弛信息。PMP 不是取代变分法,而是包含变分法作为特例。

从变分法到 Hamilton-Jacobi-Bellman

除了 PMP(沿单条轨迹的必要条件),变分法还有另一条推广路径:Hamilton-Jacobi-Bellman (HJB) 方程

Hamilton-Jacobi 方程的来源:考虑最优值函数 \(V(x,t) = \min_{y:[t,T]} J[y]\)(从 \((x,t)\) 出发的最优代价)。利用 Bellman 最优性原理(最优轨迹的任意后段仍是最优的),可以推导出 \(V\) 满足的 PDE:

\[\frac{\partial V}{\partial t} + \min_u \left[L(x,u,t) + \nabla_x V \cdot f(x,u)\right] = 0\]

这就是 HJB 方程。它是一个一阶非线性 PDE(在控制设定下),与 PMP 的 ODE 系统形成对偶。

PMP 与 HJB 的关系

\[\lambda(t) = \nabla_x V(x^*(t), t)\]

PMP 的共态 \(\lambda(t)\) 正是 HJB 值函数 \(V\) 沿最优轨迹的梯度。PMP 是 HJB 的**特征线方程**——就像一阶 PDE 的特征线方法一样。

特征 PMP HJB
方程类型 ODE(沿单条轨迹) PDE(全状态空间)
条件性质 必要条件 充分条件(如果光滑解存在)
信息量 只给出最优轨迹 给出从任意初始状态的最优策略
计算难度 \(O(n)\)(状态维度) \(O(N^n)\)(维数灾难)
适用场景 开环轨迹优化 闭环反馈控制、RL

变分法中的 Hamilton-Jacobi 理论:即使在经典变分法中(不涉及控制),也有 Hamilton-Jacobi 方程。对标准泛函 \(J = \int L\,dx\),定义"作用函数" \(S(x,y) = \min_{y:[a,x]} J[y]\)(从左端到点 \((x,y)\) 的最优代价),则 \(S\) 满足

\[\frac{\partial S}{\partial x} + H\left(x, y, \frac{\partial S}{\partial y}\right) = 0\]

其中 \(H(x,y,p) = py' - L(x,y,y')\)\(y'\) 通过 \(p = L_{y'}\) 反解)。这就是 Hamilton-Jacobi 方程——经典力学中的核心方程之一。

跨领域类比:HJB 方程在最优控制中的地位,类似于 Bellman 方程在强化学习中的地位。事实上 Bellman 方程就是 HJB 的离散时间版本。当你学习 RL 中的 Q-learning 或策略梯度时,你使用的 Bellman 递推

\[V(s) = \min_a [r(s,a) + \gamma V(s')]\]

正是 HJB 方程在离散 MDP 上的特殊化。变分法 → HJB → Bellman → RL 是一条清晰的数学传承链。


本章小结

知识点 核心内容 难度 工程应用
泛函与变分 函数空间中的"导数" 所有后续推导的基础
EL 方程 \(L_y - \frac{d}{dx}L_{y'} = 0\) 机器人动力学、轨迹优化
变分引理 \(\int f\eta = 0 \Rightarrow f \equiv 0\) EL 推导的关键步骤
Beltrami 恒等式 \(L - y'L_{y'} = C\)\(L\) 不含 \(x\) ⭐⭐ 快速降阶求解
等周约束 增广 Lagrangian \(\tilde{L} = L + \lambda G\) ⭐⭐ 约束优化的原型
横截条件 自由端点 → \(L_{y'} = 0\) ⭐⭐ PMP 横截条件的前身
Hamilton 原理 \(\delta S = 0\) → Lagrange 方程 ⭐⭐ 机器人动力学方程来源
Legendre 变换 \((q,\dot{q}) \to (q,p)\) ⭐⭐ PMP 框架的前奏
辛结构 Hamilton 方程 \(\dot{z} = J\nabla H\) ⭐⭐⭐ 辛积分器、PMP 结构
Noether 定理 连续对称性 ↔ 守恒量 ⭐⭐⭐ 辛积分器、守恒律
Legendre 条件 \(L_{y'y'} \geq 0\) ⭐⭐ 验证极小性
Jacobi 条件 无共轭点 → 弱极小 ⭐⭐⭐ 轨迹优化局部最优验证
Weierstrass E-函数 \(\mathcal{E} \geq 0\) → 强极小 ⭐⭐⭐ 理论完备性保证
Ritz/Galerkin 法 直接法:有限维近似 ⭐⭐⭐ 有限元、轨迹优化的基础
EL → PMP 桥梁 \(\partial H/\partial u = 0\)\(\max_u H\) ⭐⭐⭐ 最优控制理论核心

累积项目:本章新增模块

项目:从零构建最优控制求解器

本章新增:变分 BVP 求解模块 - 实现 EL 方程的符号推导(SymPy euler_equations) - 用 scipy.integrate.solve_bvp 数值求解最速降线 BVP - 用 Ritz 法近似求解简单变分问题并与解析解对比 - 验证解析解(摆线)与数值解的一致性 - 绘制解曲线和代价函数收敛图

# 累积项目代码骨架
import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

# === 模块 1:EL 方程符号推导 ===
from sympy import symbols, Function, diff, simplify
from sympy.calculus.euler import euler_equations

x = symbols('x')
y = Function('y')(x)
# 例:最速降线的 Lagrangian
L = ((1 + diff(y, x)**2) / (2 * 9.81 * y))**0.5
# EL 方程(符号形式)
# el_eq = euler_equations(L, y, x)  # SymPy 自动推导

# === 模块 2:BVP 数值求解 ===
def brachistochrone_ode(x, y):
    """最速降线 ODE 系统(由 EL 方程化为一阶系统)
    状态: y[0] = y(x), y[1] = y'(x)
    由 Beltrami 恒等式 y(1+y'^2) = 2r 推导 y'' 表达式
    """
    y_val, yp = y[0], y[1]
    if y_val < 1e-10:
        y_val = 1e-10  # 避免除零(x=0 处的奇异性)
    ypp = (1 + yp**2) / (2 * y_val)
    return np.array([yp, ypp])

def brachistochrone_bc(ya, yb):
    """边界条件: y(0)=0 (近似为小值), y(x1)=y1"""
    return np.array([ya[0] - 1e-6, yb[0] - 1.0])

# 初始猜测与求解
x_mesh = np.linspace(0, 1, 100)
y_init = np.zeros((2, x_mesh.size))
y_init[0] = np.linspace(1e-6, 1.0, x_mesh.size)  # y 的初始猜测
y_init[1] = 1.0  # y' 的初始猜测

sol = solve_bvp(brachistochrone_ode, brachistochrone_bc, x_mesh, y_init)

# === 模块 3:Ritz 法近似 ===
def ritz_solve_simple(N_basis=5):
    """用 Ritz 法求解 min ∫[(y')^2 + y^2 - 2xy] dx, y(0)=y(1)=0
    基函数: sin(k*pi*x), k=1,...,N
    """
    from numpy.linalg import solve as lin_solve
    # 组装刚度矩阵 K 和载荷向量 f
    K = np.zeros((N_basis, N_basis))
    f_vec = np.zeros(N_basis)
    for i in range(N_basis):
        for j in range(N_basis):
            ki, kj = (i+1)*np.pi, (j+1)*np.pi
            # ∫ sin(ki*x)*sin(kj*x) dx = 0.5*delta_{ij}
            # ∫ cos(ki*x)*cos(kj*x)*ki*kj dx = 0.5*ki*kj*delta_{ij}
            if i == j:
                K[i,j] = 0.5 * (ki**2 + 1)
            # else: K[i,j] = 0 (正交性)
        # ∫ x*sin(ki*x) dx (分部积分计算)
        ki = (i+1)*np.pi
        f_vec[i] = (-1)**(i+2) / ki  # = (-1)^{i+1}/ki
    c = lin_solve(K, f_vec)
    return c

# 后续章节将扩展此框架到 PMP 间接法和 DDP

下一章扩展方向: - 模块 4:PMP 必要条件的验证(共态方程数值积分) - 模块 5:直接配点法(Collocation)实现 - 模块 6:与 CasADi 求解器结果对比

# === 模块 4:Ritz 法可视化对比 ===
def ritz_visualization():
    """对比不同基函数数目的 Ritz 近似"""
    x = np.linspace(0, 1, 200)

    plt.figure(figsize=(10, 6))
    for N in [1, 2, 5, 10]:
        c = ritz_solve_simple(N)
        y_approx = sum(c[k] * np.sin((k+1)*np.pi*x) for k in range(N))
        plt.plot(x, y_approx, label=f'Ritz N={N}')

    plt.xlabel('x')
    plt.ylabel('y(x)')
    plt.title('Ritz 法收敛性:基函数数目的影响')
    plt.legend()
    plt.grid(True)
    plt.show()

延伸阅读

资源 难度 内容
Gelfand & Fomin《Calculus of Variations》Dover ⭐⭐ 经典入门,230页精炼
Liberzon《Calculus of Variations and Optimal Control》(免费 PDF) ⭐⭐⭐ 唯一单本贯通 CoV→PMP→HJB
Lanczos《The Variational Principles of Mechanics》Dover ⭐⭐ 物理直觉最佳
Kirk《Optimal Control Theory》Dover ⭐⭐ 工程视角入门
van Brunt《The Calculus of Variations》Springer ⭐⭐⭐ 全面含 Noether 专章
Sussmann & Willems "300 Years of Optimal Control" (IEEE 1997) ⭐⭐ 历史综述必读
3Blue1Brown "The Brachistochrone" (YouTube) 最佳视觉入门
Liberzon 讲义 PDF: liberzon.csl.illinois.edu/teaching/cvoc.pdf ⭐⭐⭐ 免费研究生教材
Marsden & West "Discrete Mechanics and Variational Integrators" (Acta Numerica, 2001) ⭐⭐⭐⭐ 离散变分力学经典
Mellinger & Kumar "Minimum Snap Trajectory Generation" (ICRA, 2011) ⭐⭐ 四旋翼轨迹规划的变分基础
Courant & Hilbert《Methods of Mathematical Physics》Vol. I, Ch. IV ⭐⭐⭐ 变分法的数学物理视角
Giaquinta & Hildebrandt《Calculus of Variations I》Springer ⭐⭐⭐⭐ 现代变分法全面参考

🔧 故障排查手册

症状 可能原因 排查步骤 相关章节
EL 方程解不满足端点条件 常数选取错误或无解 1. 检查问题是否有解(存在性定理) 2. 验证 BVP 可解性条件 3. 检查积分常数 4. 对多解问题检查所有分支 §1.2, §1.3
solve_bvp 不收敛 初始猜测太差 1. 用解析近似作初始猜测 2. 增加网格点 3. 对奇异问题做正则化(如将 \(y(0)=0\) 改为 \(y(0)=\epsilon\)) 4. 尝试 continuation 方法 §1.3, §1.8
Legendre 条件不满足 EL 解不是极小值 1. 检查是否为极大值(换符号重新求解) 2. 寻找其他极值曲线 3. 检查问题是否有极小值(如泛函无界则无极小) §1.7
Beltrami 恒等式等号右边为零 特殊情形,需要额外分析 1. 检查 \(L\) 是否关于 \(y'\) 齐次 2. 考虑退化解 3. 用完整 EL 方程验证 §1.6
约束问题乘子 \(\lambda\) 无法确定 约束退化或异常情形 1. 检查约束是否"活跃" 2. 验证正则条件(\(\delta K \neq 0\)) 3. 考虑异常乘子情形(\(\lambda_0 = 0\) §1.4
Ritz 法收敛缓慢 基函数选择不当 1. 增加基函数数目 2. 换用更适合问题结构的基函数 3. 检查解是否有奇异性(角点、边界层) §1.8
数值解在端点附近振荡 边界层效应或 Gibbs 现象 1. 在边界附近加密网格 2. 使用自适应基函数 3. 检查是否需要弱解理论 §1.8

方法间关系图:从变分法到现代最优控制

变分法(1696-1744, Bernoulli-Euler-Lagrange)
├── 无约束: Euler-Lagrange 方程 ────────────────────┐
│   ├── 力学: Hamilton 原理 → Lagrange/Hamilton 方程   │
│   ├── 几何: 测地线方程 → SLAM/规划中的流形优化       │
│   └── 直接法: Ritz/Galerkin → FEM → 轨迹优化       │
│                                                      │
├── 等周约束: Lagrange 乘子法 ─────────────┐          │
│                                          │          │
│   控制约束 u∈U 的本质困难:               │          │
│   ∂H/∂u=0 在边界上失效!                │          │
│   └── 需要新工具 ──────────────────────────┤          │
│                                          ▼          ▼
│                              Pontryagin 极大值原理 (PMP, 1956-1962)
│                              │ 必要条件, 沿单条轨迹, ODE
│                              │
│                              ├── 离散化 ──→ DDP/iLQR (1970-2012)
│                              │             └── Crocoddyl, OCS2
│                              │
│                              └── 对偶 ──→ Hamilton-Jacobi-Bellman (HJB)
│                                          │ 充分条件, 全状态空间, PDE
│                                          │
│                                          ├── 离散 ──→ Bellman 方程 / DP (1957)
│                                          │             │
│                                          │             └── 近似 ──→ RL (1990s-)
│                                          │
│                                          └── 数值 ──→ Level Set / DeepReach
├── 直接法分支:
│   ├── Ritz (1908) → FEM (1960s) → 工程仿真
│   ├── Galerkin (1915) → 谱方法 → 高精度 PDE 求解
│   └── 配点法 (Collocation) → 直接转录 → CasADi/IPOPT/Drake
└── 统一视角: PMP 是 HJB 的特征线; HJB 是 PMP 的"包络面"
    共态与值函数梯度关系: lambda(t) = grad_x V(x*(t), t) ← 两条路殊途同归

变分法的历史时间线

年代 人物 贡献 数学意义
1696 Johann Bernoulli 提出最速降线问题 变分法的诞生契机
1697 Jakob Bernoulli 通用变分方法 从特殊技巧到一般方法
1697 Newton 匿名投稿解答 "从爪子认出狮子"
1744 Euler 系统化 EL 方程 建立理论框架
1755 Lagrange \(\delta\) 记号、纯分析方法 消除几何直觉依赖
1786 Legendre 二阶条件(有缺陷) 开始研究充分条件
1834 Hamilton Hamilton 原理、正则方程 统一力学形式
1837 Jacobi 共轭点、Jacobi 方程 修补 Legendre 的不足
1870 Weierstrass E-函数、强极小 严格化充分条件
1879 du Bois-Reymond 变分基本引理推广 弱解理论的先驱
1908 Ritz Ritz 近似法 直接法的开端
1915 Galerkin 加权残差法 FEM 的理论基础
1918 Noether 对称性与守恒量 20 世纪最重要的数学定理之一
1956 Pontryagin 极大值原理 从 EL 到 PMP 的飞跃
1957 Bellman 动态规划 HJB 的离散版本
2001 Marsden & West 离散变分力学 变分积分器

这条时间线跨越三个世纪,每一步都建立在前人的基础上。学习变分法不只是学习一套技术,更是体会数学思想的演化——从具体问题到一般理论,从必要条件到充分条件,从连续到离散,从无约束到有约束。


本章常见误解汇总

误解 正确理解
EL 方程的解一定是极小值 EL 方程只是一阶必要条件(类似 \(f'(x)=0\)),可能对应极大值、极小值或鞍点;需要二阶条件(Legendre、Jacobi、Weierstrass)确认
变分法中的 \(\delta y\) 是一个无穷小量 \(\delta y = \varepsilon \eta(x)\) 是一个有限函数乘以参数 \(\varepsilon\),"第一变分为零"指的是泛函关于 \(\varepsilon\) 的导数在 \(\varepsilon=0\) 处为零
最小作用量原理说物理系统总是使作用量最小 准确表述是作用量取**驻值**(stationary),不一定是最小值;在某些情况下(如共轭点之后)实际轨迹对应的是鞍点
Beltrami 恒等式和 Noether 定理是同一回事 Beltrami 恒等式是 \(L\) 不显含自变量时的特殊第一积分,Noether 定理是任意连续对称性与守恒量之间的一般对应;前者是后者的特例
变分法只能处理固定端点问题 横截条件(transversality condition)允许处理自由端点、可变终端时间等情况,是变分法的自然推广
Ritz 法和 Galerkin 法是近似方法,不严格 在适当的函数空间和收敛性假设下,这些方法有严格的收敛性理论;有限元方法(FEM)正是基于 Galerkin 框架
EL 方程可以处理所有优化问题 EL 方程假设控制/决策变量无界且解足够光滑;当存在约束(如 \(u \in U\))时,需要 PMP;当解不光滑时,需要弱解理论
等周问题的 Lagrange 乘子 \(\lambda\) 只是数学技巧 \(\lambda\) 具有明确的物理意义:它衡量约束条件松弛一个单位时目标泛函的最优值变化量(影子价格)

向后展望:下一章预告

本章建立了变分法的完整理论体系:

  • 必要条件:EL 方程(一阶)、第二变分相关条件(二阶)
  • 充分条件:Legendre + Jacobi(弱极小)、Weierstrass(强极小)
  • 计算工具:Beltrami 恒等式(降阶)、Noether 定理(守恒量)、Ritz/Galerkin(数值)
  • 物理联系:Hamilton 原理、Lagrange 力学、辛结构

但本章有一个根本局限:控制变量没有约束

在真实机器人系统中,控制永远是有界的:\(|u| \leq u_{\max}\)(推力饱和)、\(\tau \in [\tau_{\min}, \tau_{\max}]\)(力矩限制)。当最优控制恰好在约束边界上时(如 bang-bang 控制),EL 方程的核心步骤 \(\partial H / \partial u = 0\) 失效——因为极值点不在内部!

这正是 Pontryagin 极大值原理(PMP,下一章)要解决的问题:用"在约束集 \(U\) 上极大化 Hamiltonian"替代"Hamiltonian 对 \(u\) 的梯度为零",从而处理任意紧集 \(U\) 上的控制约束——包括非凸约束、甚至离散控制集。

PMP 还解决了变分法的另一个局限:自由终端时间。在变分法中,积分区间 \([a,b]\) 是固定的,虽然本章 1.4.5 节的横截条件可以处理可变端点,但方法较为繁琐。在最优控制中,终端时间 \(T\) 本身经常是优化变量——例如"最短时间到达目标"(时间最优控制)。PMP 通过附加条件 \(H(t_f) = 0\) 简洁地确定最优 \(T\),将固定时间和自由时间问题统一到同一理论框架中。

下一章的核心公式——

\[H(x^*, \lambda^*, u^*, t) \geq H(x^*, \lambda^*, u, t), \quad \forall u \in U\]

——就是变分法 \(\partial H/\partial u = 0\) 在"\(u\) 有约束"情形下的正确推广。

学习建议:在进入下一章之前,请确保你能独立完成以下三件事: 1. 对给定 Lagrangian 写出 EL 方程(包括多维情形) 2. 识别 \(L\) 不含 \(x\) 或不含 \(y\) 的特殊情形并应用 Beltrami 恒等式或动量守恒 3. 完整求解至少一个经典变分问题(最速降线或悬链线)

如果对以上任何一项有困难,请回到对应章节重新学习后再继续。变分法是最优控制的地基——地基不牢,PMP 和 DDP 的大厦无法建立。