简介面向需要求解Hamilton-Jacobi-Bellman(HJB)方程的研究人员与工程开发者这份MATLAB工具包提供了一维空间与控制设定下的数值求解实现。核心函数HJB.m、optimization.m与reachableset.m分别负责主迭代、单步优化与可达集计算配合createsets.m、refine.m等辅助脚本支持用户在给定时间区间、空间步长与边界条件下快速完成动态规划与最优控制问题的离散求解和实验复现。资源打包为zip格式共16个文件以11个m脚本为主另含2张示意图、1份HTML与1份Markdown说明文档以及1个演示用GIF动图整体仅124KB结构精简适合直接阅读代码或在此基础上扩展。目前已有1975人学习下载适合具备一定MATLAB基础、希望理解HJB方程数值解法或开展最优控制教学实验的读者使用。1. HJB方程数值解为什么总在最后一公里卡住你Hamilton-Jacobi-Bellman 方程也就是常说的 HJB 方程是最优控制问题的连续时间表述核心。很多从业者第一次接触它不是在教科书里而是在做轨迹规划、库存控制、期权定价或者自动驾驶决策时发现自己辛辛苦苦推导的反馈控制律根本没法落地——因为解析解只在极简单的线性二次型问题里存在现实问题几乎全都要靠数值解硬啃。我见过不少团队把强化学习框架搭得风生水起最后落到连续状态空间的值函数逼近时才发现 HJB 方程的数值格式比想象中敏感得多网格密了算不动网格稀了精度崩时间步长调大了直接发散。这个标题让你进来的原因大概率就是你正卡在这条路上需要算出某个最优控制问题的值函数跑出一条能用的反馈策略。这篇文章的目的只有一个把 HJB 方程的数值解从「理论研究」掰成「工程活」。我会先讲清楚它到底在解什么再给出有限差分、半拉格朗日离散加策略迭代这套最常见的落地路径然后按我的血泪经验列出最容易翻车的几个坑。无论你是手写求解器还是准备基于某些现成库二次开发这篇文章都能帮你少走几轮弯路。2. HJB方程从最优性原理到可离散的PDE先把要解的对象写对2.1 确定性最优控制问题里的 HJB 长什么样先限定在确定性系统。考虑一个连续时间控制系统[ \dot{x}(t) f(x(t), u(t)) ]代价函数是从当前时刻 ( t ) 到终端时刻 ( T ) 的折扣或终值积分[ V(x, t) \inf_{u} \left{ \int_t^T L(x(s), u(s)) e^{-\lambda(s-t)} ds \Phi(x(T)) \right} ]其中 ( L ) 是运行代价( \lambda ) 是折扣率( \Phi ) 是终端代价。在这个设定下动态规划原理给出值函数满足的偏微分方程[ -\frac{\partial V}{\partial t} - \lambda V \min_{u \in U} \left{ L(x, u) abla_x V \cdot f(x, u) \right} 0 ]这就是 HJB 方程的一般形式。它是个一阶非线性偏微分方程而且在绝大多数问题里粘性解不光滑经典偏微分方程理论那一套很多都不能直接用。做数值解之前必须知道自己在解粘性解意义下的方程否则后面构造的格式稳定性没有任何保证。2.2 稳态 HJB最常被实际项目采用的形式如果系统在无限时域上运行代价函数不再是时间 ( t ) 的函数而是只依赖于状态 ( x )那么这个方程退化为一个不含时间导数的椭圆型 PDE。金融里常见的衍生品定价、平稳库存控制、无限时域 LQR 的推广用的都是这个稳态形式[ \lambda V(x) \min_{u \in U} \left{ L(x, u) abla V(x) \cdot f(x, u) \right} ]注意这里的 ( \lambda ) 不再是折扣率那么简单它必须大于零否则方程可能没有粘性解或者解无界。很多初次接触 HJB 的开发者在这里就栽了跟头直接把折扣率设成零结果迭代怎么算都是发散或震荡。实际落地时( \lambda ) 的取值直接影响解的稳定性和迭代收敛速度。2.3 一个能跟着思考的最小例子一维积分器为了不让你被符号淹没我拿出一个最简单但完整的例子一维状态 ( x )控制 ( u \in [-1, 1] )系统为 ( \dot{x} u )运行代价 ( L(x, u) \frac{1}{2} x^2 \frac{1}{2} u^2 )折扣率 ( \lambda 0.05 )。这是一个标准的无限时域调节问题。它的 HJB 方程为[ 0.05 V(x) \min_{u \in [-1,1]} \left{ \frac{1}{2} x^2 \frac{1}{2} u^2 V(x) \cdot u \right} ]这个方程存在精确解只是需要分段讨论控制约束是否起作用。用这个例子做数值解的出口测试非常合适后面我会反复用到它。2.4 从变分法到上风差分为什么直接离散是错的这里的关键数值问题是对 ( V(x) ) 做中心差分也就是 ( \frac{V_{i1}-V_{i-1}}{2\Delta x} )看起来相当自然但在这个 HJB 方程里大概率直接炸掉。原因在于方程里的哈密顿量 ( \min_{u} { L V(x) \cdot f } ) 对 ( V(x) ) 的依赖可能是非线性的而中心差分会引入不自洽的数值耗散导致格式不稳定。我见过一个项目用中心差分解一个带强约束的 HJB 问题残差振荡了上千次迭代都不收敛最后换成上风差分才稳定下来。你需要按控制的方向来判断信息传播的方向也就是所谓「迎风」或者「上风」思想。具体操作会在下一章的离散格式里展开。3. 把 HJB 变成可计算的离散方程有限差分格式与策略迭代骨架3.1 选择网格均匀网格还是自适应网格先明确一点HJB 数值解几乎永远在某个有界计算域 ( \Omega \subset \mathbb{R}^n ) 上进行。你不可能在无穷状态空间上铺一层均匀网格。常见做法是取一个足够大的有界盒子 ( [x_{\min}, x_{\max}]^n )超出边界时用自定义的边界条件代替真实行为。这个「足够大」需要靠实验去确认后面在避坑章节我会专门讲怎么判断网格边界是否影响了解。自适应网格例如基于梯度误差估计的网格加密可以显著减少网格点数量。但引入网格重构后值函数在不同网格层之间的插值误差也要计入总误差这在量产代码里很难控制。我的一般建议是第一版求解器先做均匀网格把格式和迭代逻辑跑通再考虑自适应加密。盲目追自适应网格会导致你分不清是格式问题还是网格问题。3.2 一阶上风差分格式怎么写把状态空间离散成节点 ( x_i )节点间距为 ( \Delta x )。对一维问题上风差分的意思是根据控制 ( u ) 和漂移函数 ( f ) 的符号决定从哪一侧取差分。在最小化框架 ( \min_{u} { L(x,u) V(x) f(x,u) } ) 中若对某个给定控制 ( u )( f(x,u) 0 )则状态流从左向右传播( V(x) ) 应该用后向差分 ( \frac{V_i - V_{i-1}}{\Delta x} )若 ( f(x,u) 0 )则用前向差分 ( \frac{V_{i1} - V_i}{\Delta x} )。把这两个情况写进最小化过程得到的离散方程是[ \lambda V_i \min_{u} \left{ L(x_i,u) \frac{V_i - V_{i-1}}{\Delta x} f^(x_i,u) \frac{V_{i1} - V_i}{\Delta x} f^-(x_i,u) \right} ]其中 ( f^ \max(f,0) )( f^- \min(f,0) )。这个格式被称为一致单调格式的一种最简单实现。它保证了信息沿特征方向传播避免中心差分带来的伪振荡。3.3 策略迭代的三步循环策略评估、策略改进、收敛判断策略迭代是解 HJB 方程最实用的不动点方法之一尤其适合控制空间离散化后的情形。它分三步走策略评估固定当前策略 ( u^{(k)} )把 HJB 方程里的min去掉求解一个线性方程。策略改进对每个状态节点根据当前值函数梯度重新计算最优控制 ( u^{(k1)}(x_i) )。收敛判断检查 ( u^{(k1)} ) 是否与 ( u^{(k)} ) 一致或者更稳健地检查值函数变化量小于阈值。这个循环的最大优点是只要初始策略别太离谱它往往比直接价值迭代收敛得快得多。代价是每轮策略评估要求解一个大规模线性方程组在二维以上状态空间里这个求解器的性能就是瓶颈。3.4 最小可运行的策略迭代Python骨架下面这份代码只有一个作用让你看到策略迭代在连续状态 HJB 里是怎么组织起来的。它针对前面的一维积分器例子用numpy实现不走任何花哨优化。import numpy as np # 参数设定 x_min, x_max -5.0, 5.0 N 201 # 网格节点数 dx (x_max - x_min) / (N - 1) lam 0.05 # 折扣率 u_max 1.0 # 控制约束 u_grid np.linspace(-u_max, u_max, 21) # 控制离散候选集 # 初始化网格与值函数 x_grid np.linspace(x_min, x_max, N) V np.zeros(N) # 值函数初始猜测 u_policy np.zeros(N) # 初始策略 # 边界外的值处理这里简单做零值外推 # 实际工程中常用线性外推或自定义边界条件 def f(x, u): return u def L(x, u): return 0.5 * x**2 0.5 * u**2 # 策略迭代主循环 for it in range(200): # 策略评估解线性方程组 A np.zeros((N, N)) b np.zeros(N) for i in range(N): u u_policy[i] fp max(f(x_grid[i], u), 0.0) fm min(f(x_grid[i], u), 0.0) # 上风差分对应的线性算子 # lambda V L (V_i - V_{i-1})/dx * fp (V_{i1} - V_i)/dx * fm # 整理成 A V b A[i, i] lam (fp - fm) / dx if i 0: A[i, i-1] -fp / dx if i N - 1: A[i, i1] fm / dx b[i] L(x_grid[i], u) V_new np.linalg.solve(A, b) # 策略改进逐点最小化 policy_changed False for i in range(N): best_u u_policy[i] best_val None for u in u_grid: fp max(f(x_grid[i], u), 0.0) fm min(f(x_grid[i], u), 0.0) # 离散导数用上风格式 if fp 0: dv (V_new[i] - V_new[i-1]) / dx if i 0 else 0.0 else: dv (V_new[i1] - V_new[i]) / dx if i N - 1 else 0.0 # 更准确的做法上下风同时贡献 dv_up (V_new[i] - V_new[i-1]) / dx if i 0 else 0.0 dv_down (V_new[i1] - V_new[i]) / dx if i N - 1 else 0.0 ham L(x_grid[i], u) dv_up * fp dv_down * fm if best_val is None or ham best_val: best_val ham best_u u if abs(best_u - u_policy[i]) 1e-12: policy_changed True u_policy[i] best_u V V_new if not policy_changed: print(f策略在第 {it} 次迭代稳定) break代码逻辑并不复杂。策略评估阶段把上风差分写进矩阵A构造出线性方程组 ( A V b )策略改进阶段在控制候选集上做枚举最小化把使哈密顿量最小的控制赋值给u_policy。需要注意dv_up与dv_down的计算方式我在这段示例里其实同时保留了两个方向的分量与纯上风思想略有差异。对一维问题这样处理没有大问题但严格一致单调格式应该只取控制流方向对应的一侧。实际使用中建议先按严格上风做如果网格不够密出现消耗振荡再调整格式设计。参数方面N201在一维问题上提供了约 0.05 的网格间距对大多数平滑问题足够u_grid取 21 个离散值更精细的控制网格会提高策略精度但增加计算量。lam0.05保证方程非奇异避免纯边值问题的求解困难。3.5 为什么策略评估要解线性方程组而不是显式迭代你可能注意到我直接用np.linalg.solve而不是价值迭代里的显式推进。原因很简单策略评估是线性的直接法在二维以下、节点数低于十万时速度可观且无条件稳定。如果你用V A_inv b之外的方式比如 Gauss-Seidel 迭代也完全可以但要注意迭代终止准则。直接法在这里最省心。4. 走向高维和复杂控制约束半拉格朗日离散与高阶修正4.1 半拉格朗日方法的基本思想把 PDE 还原成特性线问题标准的有限差分方法直接在网格上用导数近似代替方程中的导数项。半拉格朗日Semi-Lagrangian方法则换了一个角度对于确定性控制问题值函数沿最优轨迹满足动态规划原理因此可以用「从一个网格点出发沿特性线追踪到前一个时刻」的方式来构造离散方程。在数值实现里这通常体现为二维或三维空间中的插值回退。对稳态 HJB 方程半拉格朗日离散的一般形式可以写为[ V(x) \min_{u} \left{ e^{-\lambda \Delta t} V(x f(x,u) \Delta t) \int_0^{\Delta t} e^{-\lambda s} L(x(s), u) ds \right} ]然后你把 ( x f(x,u) \Delta t ) 落到相邻网格之间用线性插值取 ( V ) 在这个非网格点处的近似值。这种方法的天然优势是即使在很粗的网格上它对特征方向的追踪也比纯上风差分更准因为它在做「沿方向走一步再插值」而不是生硬地区分前向和后向差分。4.2 半拉格朗日离散怎么写成代码在实现时半拉格朗日格式反而比有限差分更直观对每个节点枚举候选控制计算一步回退位置然后在网格上做插值。下面这段是核心循环框架我对前面的一维例子做了替换。import numpy as np # 参数 x_min, x_max -5.0, 5.0 N 201 dx (x_max - x_min) / (N - 1) lam 0.05 dt 0.8 * dx / u_max # 满足 CFL 条件的回退步长 u_max 1.0 u_grid np.linspace(-u_max, u_max, 21) x_grid np.linspace(x_min, x_max, N) def L(x, u): return 0.5 * x**2 0.5 * u**2 def f(x, u): return u # 线性插值函数超出边界时用边界节点的值外推 def interp_V(V, xq): idx (xq - x_min) / dx i0 int(np.floor(idx)) if i0 0 or i0 N - 1: return V[0] if i0 0 else V[-1] ratio idx - i0 return (1 - ratio) * V[i0] ratio * V[i0 1] # 策略迭代 V np.zeros(N) u_policy np.zeros(N) for it in range(300): V_prev V.copy() # 策略评估因为 V 出现在右侧插值项里用迭代法求解隐式方程 for _ in range(50): V_new np.zeros(N) for i in range(N): u u_policy[i] x_next x_grid[i] f(x_grid[i], u) * dt V_new[i] np.exp(-lam * dt) * interp_V(V, x_next) dt * L(x_grid[i], u) V V_new # 策略改进 policy_changed False for i in range(N): best_u u_policy[i] best_val None for u in u_grid: x_next x_grid[i] f(x_grid[i], u) * dt val np.exp(-lam * dt) * interp_V(V, x_next) dt * L(x_grid[i], u) if best_val is None or val best_val: best_val val best_u u if abs(best_u - u_policy[i]) 1e-12: policy_changed True u_policy[i] best_u if not policy_changed: print(f半拉格朗日策略迭代在第 {it} 次稳定) break这段代码里最关键的是interp_V半拉格朗日格式能不能保持单调和收敛全看这一步插值。用线性插值是最稳妥的选择因为高阶插值如二次样条在某些情况下会破坏单调性导致数值解偏离真解。如果网格够细高阶插值能换来精度收益但你必须同时验证单调性没有被打坏。dt的选择要满足 CFL 条件常见做法是让 ( \Delta t \le \Delta x / \max |f| )。我在这里取0.8 * dx / u_max留了一点余量。如果dt太大会导致回退一步跨过多个网格线性插值误差急剧上升最终解出现波浪状跳动。4.3 控制约束和非光滑哈密顿量的处理带有控制约束的 HJB 方程哈密顿量经常不可微min算子造成 kink。数值格式必须能处理 kink否则策略改进步会一直在两个候选控制造振荡。一个经验法则是把控制空间离散得足够密但不追求连续优化。离散控制空间与连续控制空间之间的误差是 ( O(\Delta u^2) ) 量级通常可以接受如果要更精确可以在策略改进步用二次插值逼近最优控制而不是单纯在离散点里选。4.4 为什么半拉格朗日方法与强化学习的连接是你该关注的许多强化学习算法在连续状态空间中用函数逼近或查表法估计价值函数HJB 方程的半拉格朗日离散在结构上与这类方法天然接近。如果你后续想比较 HJB 数值解和近似动态规划结果半拉格朗日解法是一个良好的中间参照系它能用简单网格适配高维问题的特征方向而有限差分格式在高维下的构造复杂度迅速上升。5. 数值解避坑边界、初值、迭代发散与网格密度的五个实操问题5.1 计算域边界导致解被「粘住」现象值函数在边界附近出现明显的不自然平直或大幅跳跃而且把网格往外扩后内部解也变了。原因边界条件选得不合理。最常见的错误是把边界节点上的值固定为零或者采用零导数外推两者都会产生边界层。解决对稳态 HJB 问题最稳妥的边界设定是「把边界节点上的贴现方程单独处理」例如直接用单侧差分且假设边界外的值线性外推。在我的工程实践里推荐做法是先把计算域取到比感兴趣区域大 20% 以上再观察解靠近边界 10% 范围内是否发生异常若异常就继续扩大。判断标准很简单如果边界移动了但内部解变化小于工程允许误差边界就足够远。5.2 价值迭代或策略评估迭代不收敛现象迭代几千步后残差不再下降在一个范围内振荡或者数值疯狂增长到无穷。原因可能是贴现率太小、时间步长太大、格式不是单调的也可能是策略评估里存在周期循环策略在 A 和 B 之间来回跳值函数随之震荡。解决先查lam是否大于 0如果接近零矩阵接近奇异直接法也会不稳定。再把dt缩小一半重跑若稳定下来说明 CFL 条件被突破。如果还是振荡检查策略改进的容差是不是设得过严把1e-12放宽到1e-6常常能帮助脱离死循环。这里有个比较隐蔽的因素当最优控制不在候选控制网格点上时策略会在两个相邻候选之间抖动。解决办法是加密控制网格或者用相邻候选控制插值出更精确的最优控制。5.3 网格密度改了结果却不收敛到同一个解现象N100和N200算出来的值函数在关键区域相差 5% 以上而且N继续加大差值没有缩小趋势。原因说明离散格式本身不是一致的或者边界影响没有消除。解决先做网格收敛性测试取 ( N 200, 400, 800 ) 对比关键状态点的值。如果误差随网格加密线性下降说明一阶格式没问题如果误差不变甚至上升检查格式是否真的单调。另一个常见做法是和已知解析解对比。我见过有人把上风方向写反了网格加密后误差反而增大因为格式根本不是一致收敛的。5.4 值函数迭代初值太差导致死循环现象策略评估阶段在线性方程组里出现负对角元或者极大数值迭代直接失败。原因初值 ( V 0 ) 在大多数情况下完全合理但如果你把初值设置成常数大数比如 100而真实值函数量级是 1那么策略改进步在初期会找出一堆极端控制把迭代带进沟里。解决用零初值一般不会有问题如果确实有先验知识设置初值的量级要贴近真实值。另外一个保险做法是先跑几十轮价值迭代让值函数处于合理量级再切换到策略迭代。5.5 半拉格朗日格式在强非线性系统里出现插值耗散现象用半拉格朗日格式解一个控制项系数为 ( x^2 ) 的非线性系统网格较粗时解被「抹平」峰值明显偏低。原因线性插值本身带数值耗散当特性线路径弯曲或状态变化剧烈时这种耗散会被放大。解决加密网格可以缓解但成本高。另一个更高效的思路是把局部网格做加密或者改用二次插值但必须做单调性检验。我自己的习惯是先用粗网格线性插值算一遍找出解的高梯度区域再在这些区域做局部加密配合二次插值。这样能把精度提上去又不会引入整体计算量爆炸。6. 从 HJB 到工程验证三个一致性检查让求解器可信拿到一组看起来合理的数值解千万不要直接拿去用。数值误差在偏微分方程求解里常常伪装成「看起来合理」的结果。我每次完成 HJB 求解器后都会做下面三个检查这里把它们分享给你。第一个检查是解析解对比。对于一维积分器这类带解析解或半解析解的问题直接把数值解和精确值拉一条曲线对比计算最大绝对误差和 L2 误差。如果最大误差出现在边界多半是边界处理的问题如果出现在解的高梯度区域那就是网格分辨率不足。做完这个检查你就知道当前网格能达到的精度上限后续换到没有解析解的问题时心里有底。第二个检查是策略验证。把求出的最优策略施加到原始动力学方程上做正向仿真然后计算所得轨迹的实际代价与 HJB 值函数预测的代价对比。这两者之间的差距包含两层误解风险一是 HJB 求解误差二是仿真离散误差。如果差距小于你的工程容忍度说明求解器输出是真的可用。我常把这类验证叫做「闭环验证」因为它比单纯看残差更能说明问题。下面这段示例是验证流程的骨架# 正向闭环仿真验证 def simulate_closed_loop(x0, policy_func, T, dt_sim): x x0 cost 0.0 t 0.0 while t T: u policy_func(x) # 这里使用简单的欧拉积分实际工程可以用更高阶的 RK45 等 cost (0.5*x**2 0.5*u**2) * dt_sim * np.exp(-lam * t) x u * dt_sim t dt_sim return cost # 在多个初始状态上对比 x0_list np.linspace(-4.0, 4.0, 9) for x0 in x0_list: sim_cost simulate_closed_loop(x0, lambda x: interp_policy(x), 20.0, 0.01) hjb_cost interp_V(V, x0) print(f初始状态 {x0:.2f} | HJB 预测 {hjb_cost:.4f} | 仿真实际 {sim_cost:.4f})第三个检查是策略对参数扰动的敏感性。把运行代价里的权重系数从 1.0 改成 0.8 或者 1.2重新求解 HJB观察策略结构是否发生剧烈跳变。如果权重只改变 20%策略却完全翻转这说明你的解可能处在两个局部最优之间数值求解器在其中的行为很可疑。反过来如果策略随着权重平滑变化解就算得上稳固。这三个检查做完再上真实应用就是有底气的。在工业项目里值函数几乎是整个系统的黑匣子你给它的信任度必须经过验证来充值。我自己的习惯是任何一个 HJB 求解器进了项目必须先留出两个星期的验证期专门跑这些一致性检查而不是直接接入控制器。你省掉的这两周时间最终会以三倍的代价还回来。希望这篇笔记里给出的格式、参数和检查路径能帮你把 HJB 数值解从实验阶段推到可用状态。数值偏微分方程的路子没有银弹但每一步都验证、每个参数都问为什么这种老式做法在 HJB 这类问题上恰好是最省时间的。本文还有配套的精品资源点击获取