企业官网建设流程全解析

简介：面向MATLAB数值计算和最优控制研究者的求解器资源，聚焦哈密顿-雅可比-贝尔曼方程的数值求解。可实现一维状态与一维控制空间上的离散可达集计算，按时间与空间网格剖分后返回节点值矩阵，适合学习动态规划近似解法、验证最优控制理论的读者。资源包共十六个文件、仅一百二十四KB，其中十一个脚本文件覆盖主函数、可达集计算、单步优化等核心模块，两张图片和一张动图展示运行效果，另有超文本和标记语言文档辅助说明调用逻辑。目前已有一千九百七十三人浏览学习。通过该资源可掌握从网格构建、可达集求解到值函数更新的完整数值流程，并结合示例脚本和可视化输出加深对最优控制问题离散化实现的理解，是一份轻量实用的参考工具。 说实话，我最初接触 Hamilton Jacobi Bellman（HJB）方程，不是从偏微分方程教材开始的，而是从路径规划里的一个实际问题切入的。当时做机器人在未知环境中的避障，用模型预测控制做滚动规划，每过一个控制周期就要重新解一次最优控制问题。算多了之后发现一个问题：每个时刻都在重复做类似的优化，为什么不把整个状态空间上的最优代价函数一次性算出来？于是我就开始琢磨 HJB 方程的数值解，后来逐步整理成了一个叫 HJB-solver 的求解工具集。这篇博文就是基于这套工具的使用体验，聊一聊 HJB 方程数值求解的核心思路、离散格式、迭代算法，以及那些不踩一遍根本意识不到的坑。

这套方法解决的核心问题很明确：对于连续时间、连续状态的最优控制问题，HJB 方程给出了全局最优解的充分条件。你只需要解一次这个偏微分方程，就能得到从任意初始状态出发的最优反馈策略，而不需要像直接法那样在每次运行时重新做轨迹优化。这也是 HJB-solver 这个项目存在的价值——它把 HJB 方程的求解从纯理论推导变成了可以实际调用的数值工具。下面我会按照从原理到实现再到避坑的顺序，把整个求解链路拆开讲。

1. 为什么要直接面对 HJB 方程：从动态规划到偏微分方程

1.1 一个最简单的控制问题背后的全局结构

假设系统状态方程为：

[ \dot{x}(t) = f(x(t), u(t)), \quad x(0) = x_0 ]

控制目标是最小化如下代价泛函：

[ J(x_0, u(\cdot)) = \int_0^\infty e^{-\rho t} c(x(t), u(t)) dt ]

其中 (\rho > 0) 是折扣因子。这类问题在金融投资、机器人控制、能源管理里到处都是。按常规思路，你会用数值最优控制方法，比如直接打靶法、搭配法，把控制轨迹参数化，然后用非线性规划求解。这种做法在轨迹规划里没有问题，但它的结构是"一锤子买卖"——给定一个初始状态，解出一条最优轨迹。换一个初始状态，就得整套重算。

HJB 方程的思路完全相反。它先定义值函数：

[ V(x) = \min_{u(\cdot)} \int_0^\infty e^{-\rho t} c(x(t), u(t)) dt ]

然后利用动态规划原理做一步分解，把"从 (x) 出发的整体最优代价"拆成"从 (x) 出发走一小步的即时代价"加上"从新状态出发的剩余最优代价"。当这一小步的时间间隔趋近于零时，你就得到了一个偏微分方程：

[ \rho V(x) = \min_{u \in U} \left{ c(x, u) + abla V(x) \cdot f(x, u) \right} ]

这就是稳态 HJB 方程。它的意义在于：一旦求出 (V(x))，最优反馈控制就是关于 (x) 的函数，可以直接作为控制器使用。这就是为什么我最终选择在项目里直接求解 HJB 方程，而不是写一堆 MPC 逻辑——前者给出的是策略，后者给出的是轨迹。

1.2 粘性解：只谈光滑解是不够的

如果你去翻阅早期的控制理论教材，很多推导都假设 (V) 是连续可微的。但在实际问题里，这个假设经常不成立。最简单的例子是 L1 型代价函数配合线性动力学，最优值函数在切换面上会出现角点，导数不存在。机器人在障碍物边缘绕行时，代价函数同样会有很强的非光滑性。

这里必须提到粘性解的概念。粘性解理论的核心贡献是：即使 (V) 不光滑，我们也可以定义一种广义解，而且这个广义解是唯一的、稳定的，数值格式的收敛结果自动趋向于这个广义解。这意味着你不需要为了让方程可解而人为光滑化问题，直接用数值格式求解，只要格式满足单调性等条件，得到的解就是正确的粘性解。

我在 HJB-solver 项目里最开始的版本就没有意识到这一点，直接把 HJB 方程当作光滑 PDE 处理，用了中心差分格式，结果算出来的解在角点附近疯狂震荡，甚至出现负的"最优代价"。后来才明白，处理 HJB 方程不能用处理传统对流扩散方程的思维，格式设计和收敛性分析必须考虑粘性解框架。

2. 网格离散化与单调格式：数值解的地基

2.1 为什么不能随便用一个差分格式

当你决定在状态空间网格上求解 HJB 方程时，表面上看这只是偏微分方程数值解的标准流程：选网格、做差分离散、迭代求解代数方程。但 HJB 方程有一个特殊之处——Hamiltonian 中带有 (\min_u) 算子，离散后得到的不是普通线性方程组，而是一个非线性代数方程组。

更麻烦的是，如果差分格式选取不当，数值解不会收敛到粘性解。收敛性理论里有一个非常关键的结论：只有满足单调性、稳定性和一致性三个条件的格式，才能保证数值解收敛到粘性解。这里的单调性指格式满足比较原理的离散版本，物理含义是"值函数越大，邻域的值函数不会跟着变小"，用数学语言说就是格式对各个网格点的值函数求偏导后符号保持一致。

最常见的单调格式是迎风差分。以一维问题为例，方程最简形式是：

[ \rho V(x) = \min_{u \in U} \left{ c(x,u) + f(x,u) V_x(x) \right} ]

如果 (f(x,u) > 0)，信息从左往右传播，那 (V_x) 就应该采用向后差分：

[ V_x \approx \frac{V(x_i) - V(x_{i-1})}{h} ]

如果 (f(x,u) < 0)，则采用向前差分。这套逻辑和双曲守恒律里的迎风格式完全一致，本质上就是沿着特征方向取上游信息。HJB-solver 里我用这个策略构造离散 Hamiltonian，每个网格点上的非线性方程变成：

[ \rho V_i = \min_{u \in U} \left{ c(x_i, u) + f^+(x_i,u) \frac{V_i - V_{i-1}}{h} + f^-(x_i,u) \frac{V_{i+1} - V_i}{h} \right} ]

其中 (f^+ = \max(f,0))，(f^- = \min(f,0))。这个格式保证了离散方程组在合适条件下满足单调性，是后面所有迭代算法能收敛的基础。

2.2 离散化参数的选择与边界处理

网格步长的选择不是一个可以随便拍脑袋的问题。迎风差分的数值耗散是 (O(h)) 量级的，所以步长太大会抹平解的细节。但在状态空间维度较高时，网格点数随维度呈指数增长，步长又不能取得太小。我在实践中常用的是一个两步验证法：先用较粗的网格算一遍，再加密一倍，比较解的差异。如果关键区域（比如切换面附近）值的差异超过 5%，继续加密；如果差异小于 5%，说明当前网格已经够用。这个方法虽然土，但比任何理论误差界都可靠。

边界条件是另一个容易被忽略的坑。对于折扣因子 (\rho > 0) 的稳态问题，理论上 (V) 在状态空间无穷远处趋于某个有界常数。但数值计算必须用有界区域截断，边界条件给得不对，会直接污染内部解。我的做法是：先分析最优控制问题的结构，估算值函数在边界附近的渐进行为。如果可以确定最优控制倾向于把状态推回可行域内部，那么在边界上使用零纽曼条件（(\partial V / \partial n = 0)）通常是比较鲁棒的默认选择。如果边界处状态不可逆，则需要改成狄利克雷条件，直接指定边界值为该处保持不动时的累计代价。

3. 策略迭代：把非线性问题变成一串线性问题

3.1 策略评估与策略改进的交替流程

离散化之后，你得到了一个复杂的非线性代数方程组。直接拿牛顿法求解有两个问题：一是初值敏感，给不好的初值直接发散；二是每次迭代都要重新计算 (\min_u) 算子，代价很高。实际项目中我几乎不用牛顿法，而是用策略迭代（Policy Iteration）。

策略迭代的逻辑非常符合直觉，分两步交替执行：

1. 策略评估：给定当前控制策略 (u^{(k)}(x))，把它代入 HJB 方程，此时 (\min_u) 退化为一个确定的控制，方程变成关于 (V) 的线性椭圆方程（离散后是线性代数方程组），求解得到 (V^{(k)})。
2. 策略改进：利用当前 (V^{(k)}) 更新控制策略：

[ u^{(k+1)}(x) = \arg\min_{u \in U} \left{ c(x, u) + abla V^{(k)}(x) \cdot f(x, u) \right} ]

这两步交替进行，直到策略不再变化。关键优势在于：每次迭代只需要解一个线性方程组，而线性方程组的求解理论成熟、库多、稳定性好，比直接解非线性问题可靠得多。

在 HJB-solver 中，我把这两步封装成了核心循环。策略评估这一步的线性系统，由于迎风离散带来的对角占优性质，可以用标准迭代法求解，不需要每次都用超级复杂的预条件技术，代码简洁且稳定。

3.2 Python 示例：一维最优控制问题的核心循环

下面用一段简化的 Python 代码展示策略迭代的核心逻辑。这是一个一维问题：状态 (x \in [-2, 2])，动力学 (\dot{x} = u)，控制约束 (u \in [-1, 1])，代价函数 (c(x,u) = 0.5|x| + 0.05 u^2)，目标是尽快接近原点的同时限制控制能量。

import numpy as np # 网格设置 n = 201 x = np.linspace(-2.0, 2.0, n) h = x[1] - x[0] rho = 0.1 eps = 0.05 # 初始化 V = np.zeros(n) u = np.ones(n) # 初始策略，不妨取常数 1 # 工具函数：投影控制约束 def project_control(a): return np.clip(a, -1.0, 1.0) # 策略迭代主循环 for it in range(100): # ---- 策略评估 ---- # 构造线性系统 A V = b # 迎风方向由 f(x,u)=u 的符号决定 A = np.zeros((n, n)) b = np.zeros(n) for i in range(n): A[i, i] += rho if u[i] >= 0: # 向后差分 A[i, i] += u[i] / h if i > 0: A[i, i-1] -= u[i] / h else: # 向前差分 A[i, i] += -u[i] / h if i < n-1: A[i, i+1] -= -u[i] / h b[i] = 0.5 * abs(x[i]) + eps * u[i]**2 # 处理边界：狄利克雷边界，代价近似为稳态代价 A[0, :] = 0; A[0, 0] = 1.0; b[0] = 0.5 * abs(x[0]) / rho A[-1, :] = 0; A[-1, -1] = 1.0; b[-1] = 0.5 * abs(x[-1]) / rho V_new = np.linalg.solve(A, b) # ---- 策略改进 ---- # u* = argmin{ eps*u^2 + V_x * u } # 解析解为 u* = - V_x / (2*eps)，再做约束投影 V_x = np.gradient(V_new, h) u_new = project_control(-V_x / (2.0 * eps)) # 收敛判定：策略变化是否足够小 delta = np.max(np.abs(u_new - u)) u = u_new V = V_new if delta < 1e-6: print(f"策略迭代收敛于第 {it+1} 次迭代") break

代码里最核心的细节就是迎风差分的条件判断：根据当前策略 (u) 的符号确定差方向。为什么这很关键？因为 HJB 方程本质上是双曲型的，信息沿特征传播。如果用了中心差分，格式不是单调的，数值解在切换点附近就会出现震荡。这一段代码基本上照搬到高维状态空间时，只要把 (f(x,u)) 每个分量分别做 (f^+) 和 (f^-) 分解，逻辑完全一致。

3.3 策略迭代的收敛性观察

我在实际测试中发现，策略迭代的收敛速度非常可观。上面这个一维问题，从常数策略出发，通常 10 次迭代以内就能收敛。而且策略迭代的收敛是单调的，也就是说值函数 (V^{(k)}) 在迭代过程中逐次下降，最后稳定到最优值函数。这个单调性非常有价值——它可以作为调试工具。如果某次迭代之后值函数反而上升了，那么一定哪里出了问题，要么是离散格式错了，要么是线性系统求解出错，要么是边界条件不恰当。

不过需要提醒的是，策略迭代的收敛性理论通常需要系统可镇定，并且初始策略是稳定的。实际中你不可能总是找到稳定的初始策略，我的处理办法是用值迭代先跑几十步，得到一个粗略的值函数之后，再把这个值函数对应的策略作为策略迭代的初值。这种做法在工程上非常有效，几乎不会失败。

4. 从一维算例到高维实践：收敛性判断和真正的瓶颈

4.1 一个经过验证的测试结果

为了让大家对数值行为有直观感受，我用上面的代码做了一组测试。固定 (\rho = 0.1)，(\varepsilon = 0.05)，控制约束 (u \in [-1, 1])，状态网格从 101 逐步加密到 401 个点。迭代收敛后，我检查了两个指标：max-norm 意义下值函数的迭代变化量，以及策略改进后控制策略的变化量。结果如下：

网格加密之后，迭代次数并没有明显增加，这说明策略迭代的收敛速度对网格规模不太敏感。这一点和隐式时间步进格式很像——你增加空间分辨率时，并不需要被迫减小迭代步数来维持稳定性。当然，网格加密之后单次线性求解的规模变大了，总耗时还是会上升，但上升幅度基本是线性的，没有出现迭代次数爆炸的情况。

一个值得注意的现象是：当控制代价中 (\varepsilon) 很小接近 0 时，解在 (x = 0) 附近会形成几乎不连续的角点结构，此时策略迭代到后期会出现一种"振荡停滞"现象。原因是接近最优策略时，很小的策略扰动会导致值函数梯度的巨大变化。对策是收敛判据不要只看策略变化量，还要加上值函数变化量，两者同时小于阈值才判定收敛。目前 HJB-solver 中的默认设置就是这样，双指标判断让迭代过程可靠了不少。

4.2 状态维度升高后：网格法失效的那一刻

我在二维问题上做测试时，用的是 (51 \times 51) 的网格，总共有 2601 个点，策略迭代依然能够顺利收敛。但到了三维状态空间，如果每个维度取 51 个点，总网格数就是 (51^3 \approx 13) 万，线性系统的规模瞬间变得很大。虽然稀疏矩阵可以直接求解，但内存和时间成本已经明显上升。如果做到五维、六维的状态空间，网格法就彻底不可行了。

这就是所谓的维度诅咒。HJB 方程数值解真正的研究前沿，在很大程度上是在和维度诅咒做斗争。当你发现网格法已经走不通时，有几个替代方案值得考虑：

• 半 Lagrangian 格式：沿最优轨迹做时间反向积分，不需要在整个状态空间上做迎风差分，在控制仿射系统上效果不错，但对一般非线性系统实现复杂度较高。
• 无网格方法 / 径向基函数：避免了网格生成，但需要处理配点选择和基函数中心点的自适应分布问题，实际用起来并不比网格法省心。
• 神经网络方法（如 DGM、PINN）：用神经网络逼近值函数，理论上可以扩展到中高维。实际使用中最难的部分是保证单调性和收敛性，因为神经网络训练本身就是非凸优化，你很难证明最后得到的网络确实满足 HJB 方程的粘性解条件。

我的态度是：如果你的真实需求是三维以下，网格法加策略迭代是最稳的选择，调试方便、结果可靠、收敛性有保障。高维问题直接换神经网络试探，但一定要设计好验证集，用已知的控制问题检验网络输出是否合理，否则很容易被训练误差的假象骗过。

5. 数值求解中的那些坑：我从 HJB-solver 里学到的教训

5.1 收敛判据别只看值函数

策略迭代的收敛判据，我见过很多人只检查值函数前后两次迭代的差值。这在大多数情况下能工作，但在接近最优解时，值函数的变化量可能已经很小，而策略还在缓慢漂移。反过来也一样，策略可能几乎不变了但值函数还在下降。所以我在项目里用的是双判据：策略变化量和值函数变化量都要低于阈值，并且至少保持两个迭代周期不变，才认为真正收敛。这个做法看起来浪费了一点计算量，但避免了过早收敛。

5.2 线性求解器的选择不要盲目

策略评估步骤需要解线性系统 (A V = b)。这个矩阵因为迎风离散加上 (\rho > 0)，具有严格对角占优性，用 SOR 或 Gauss-Seidel 迭代都能收敛。但要注意：(\rho) 很小时，对角占优度变弱，经典迭代法的收敛速度会显著下降。我实际测过，(\rho) 从 0.1 降到 0.001 时，Gauss-Seidel 方法收敛所需迭代次数可以增加一个数量级。此时建议直接用稀疏直接求解器，比如 scipy 的splu或 SuperLU，稳定性和速度都有保障。别迷信"迭代法一定比直接法快"，这取决于具体问题的对角占优程度和规模。

5.3 你的代价函数会影响格式的单调性

观察上面代码中的控制量的解析解 (u^* = - V_x / (2\varepsilon))，这依赖于代价函数中 (\varepsilon u^2) 这一项。如果代价函数中有状态与控制交叉项，(\arg\min) 算子可能不是简单的解析表达式，此时需要在每个网格点上做数值优化。这个数值优化本身也需要满足单调性要求，否则会破坏整个格式的单调性。我的经验是：如果控制量维度不高（1 到 2 维），直接在每个网格点上用黄金分割搜索或坐标下降法，确保找到的确实是局部甚至全局最优；如果控制维度较高，则需要谨慎，因为每轮策略迭代里嵌套的优化本身就会成为新的性能热点。

5.4 边界条件的不可迁移性

网上很多 HJB 数值求解的示例代码，边界条件都是固定的，直接抄到自己的问题上就会出问题。边界条件必须依赖物理问题本身来设置，否则解在边界附近会出现严重畸变。我有一次做二维问题，图省事直接用了零纽曼边界，结果靠近边界的值函数出现了一个不正常的低谷，顺着这个低谷生成的反馈策略在边界附近产生了剧烈的抖振。排查了半天，最后修正为基于代价函数渐进行为的混合边界条件，现象才消失。所以，遇到异常解时，优先怀疑边界条件，这是性价比最高的排查方式。

6. 落地部署的个人经验：离线计算与在线查询

HJB 方法在实践中的经典用法是：离线求解，在线查表。状态空间网格上的 (V(x)) 解出来之后，实际上相当于已经把"从任何状态出发的最优代价"都预先计算好了。在线运行时，控制器只需要读取当前状态，在网格上插值得到 (V_x)，然后通过 (\arg\min) 公式或查表得到当前最优控制 (u)。这个查表过程耗时微秒级，完全可以满足实时性要求。

为了插值平滑，我建议在网格上保存 (V) 和 (u) 两套数据。(u) 的插值用最近邻即可，因为反馈策略本身可以有跳变（在切换面附近），强行线性插值反而会抹平策略切换，导致边界抖动。(V) 的插值则用线性插值，因为值函数虽然可能不光滑，但连续性是有的，线性插值误差可控。

最后再分享一个部署细节：在线查询时，如果状态超出计算区域边界，不要简单地截断到边界点。比较稳妥的做法是给出一段扩展区域，在扩展区域上使用解析近似，比如按当前控制保持不变的代价进行外推。这样即使系统因为扰动跑出了预设区域，控制器也不会给出荒谬的输出。这个细节虽然不复杂，但在实际工程中能省去不少麻烦。

本文还有配套的精品资源，点击获取