简介面向MATLAB数值计算和最优控制研究者的求解器资源聚焦哈密顿-雅可比-贝尔曼方程的数值求解。可实现一维状态与一维控制空间上的离散可达集计算按时间与空间网格剖分后返回节点值矩阵适合学习动态规划近似解法、验证最优控制理论的读者。资源包共十六个文件、仅一百二十四KB其中十一个脚本文件覆盖主函数、可达集计算、单步优化等核心模块两张图片和一张动图展示运行效果另有超文本和标记语言文档辅助说明调用逻辑。目前已有一千九百七十三人浏览学习。通过该资源可掌握从网格构建、可达集求解到值函数更新的完整数值流程并结合示例脚本和可视化输出加深对最优控制问题离散化实现的理解是一份轻量实用的参考工具。 说实话我最初接触 Hamilton Jacobi BellmanHJB方程不是从偏微分方程教材开始的而是从路径规划里的一个实际问题切入的。当时做机器人在未知环境中的避障用模型预测控制做滚动规划每过一个控制周期就要重新解一次最优控制问题。算多了之后发现一个问题每个时刻都在重复做类似的优化为什么不把整个状态空间上的最优代价函数一次性算出来于是我就开始琢磨 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_{i1} - 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。策略迭代的逻辑非常符合直觉分两步交替执行策略评估给定当前控制策略 (u^{(k)}(x))把它代入 HJB 方程此时 (\min_u) 退化为一个确定的控制方程变成关于 (V) 的线性椭圆方程离散后是线性代数方程组求解得到 (V^{(k)})。策略改进利用当前 (V^{(k)}) 更新控制策略[ u^{(k1)}(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, i1] - -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策略迭代收敛于第 {it1} 次迭代) 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 意义下值函数的迭代变化量以及策略改进后控制策略的变化量。结果如下网格数策略迭代次数值函数迭代末次变化量控制策略变化量101124.2e-79.1e-7201113.8e-76.5e-7401125.1e-78.2e-7网格加密之后迭代次数并没有明显增加这说明策略迭代的收敛速度对网格规模不太敏感。这一点和隐式时间步进格式很像——你增加空间分辨率时并不需要被迫减小迭代步数来维持稳定性。当然网格加密之后单次线性求解的规模变大了总耗时还是会上升但上升幅度基本是线性的没有出现迭代次数爆炸的情况。一个值得注意的现象是当控制代价中 (\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) 的插值则用线性插值因为值函数虽然可能不光滑但连续性是有的线性插值误差可控。最后再分享一个部署细节在线查询时如果状态超出计算区域边界不要简单地截断到边界点。比较稳妥的做法是给出一段扩展区域在扩展区域上使用解析近似比如按当前控制保持不变的代价进行外推。这样即使系统因为扰动跑出了预设区域控制器也不会给出荒谬的输出。这个细节虽然不复杂但在实际工程中能省去不少麻烦。本文还有配套的精品资源点击获取