微分博弈数值求解与HJI方程：追逃场景开源库实战

简介：这是一份面向博弈论、控制理论及计算数学研究者的开源微分博弈项目。项目以哈密顿-雅可比-贝尔曼-伊萨克斯（HJBI）方程为核心，展示如何用数值方法求解动态博弈中的最优策略与纳什均衡，适合研究生、算法工程师及对自动驾驶、经济建模等应用场景感兴趣的开发者学习参考。资源共含9个文件，以C语言源码为主，覆盖求解算法实现、图形可视化模块与文献管理工具，另有PDF文档说明优化方法与实验结果，并提供build.sh脚本帮助用户在Linux环境下快速编译运行。压缩包仅206KB，代码结构紧凑，便于直接阅读与二次开发。目前已有237人学习下载，可作为一个轻量而完整的入门样例，帮助读者理解HJBI方程的建模思路、数值求解流程及动态博弈结果的图形化展示方法。

1. 微分博弈的求解门槛：difgames 这个开源库替你做了哪件事

做机器人对抗、无人机追逃方向的工程师，早晚会撞上一堵墙：博弈逻辑好写，但 Hamilton-Jacobi-Isaacs（HJI）方程数值解不出来。这个方程里带 min/max 的非光滑 Hamiltonian，常规有限差分一上就振荡，很多人卡在这一步就换题了。difgames 是一个开源的微分博弈数值库，C++ 写的 HJI 求解内核，外面包了 Python 和 Julia 接口，把网格离散、时间推进、目标集处理整套都封装好。你拿到源码后不用从零推公式，直接配置场景就能跑出捕获区域、价值函数地图和最优轨迹。适合两类人：一类是做对抗控制、路径规划，需要一张可靠的捕获区地图；另一类是写论文做毕设，需要一个能改参数、能复现结果的开源基准。

2. 先把 HJI 方程啃下来：迎风差分、Lax-Friedrichs 与 CFL 这三个关键点

很多教程喜欢直接甩公式，但你不理解数值格式，后面调参翻车了连原因都猜不到。difgames 这类库的求解器跑的是 J 方程，数值内核没有多玄乎，但离散格式选错，出来的价值函数就是一片噪点。这一章我把方程长什么样、库里怎么离散、时间步长怎么定拆开讲。

2.1 先看清要解的方程：时变 HJI 与稳态变分不等式

追逃博弈里，设玩家 A 控制追逐者、玩家 B 控制逃逸者，双方都知道对方在最优对抗。定义价值函数 V(x)，表示从状态 x 出发，双方都采取最优策略时，追逐者抓到逃逸者需要的时间。这个 V 满足的是时变 HJI 方程：

∂V/∂t + min_{u∈U} max_{d∈D} { ∇V · f(x, u, d) + g(x, u, d) } = 0

其中 u 和 d 分别是追方和逃方的控制输入，f 是系统动力学，g 是累计代价。如果只关心稳态捕获时间，方程退化成 min/max 状态下的变分不等式。difgames 的经典求解对象是这个稳态解：在目标集合上 V = 0，其余区域 V 满足 Hamilton-Jacobi 条件。

难点在于 min 和 max 套在同一个 Hamiltonian 里，导致 H(∇V) 不是光滑函数。普通中心差分处理不了这种非光滑结构，强行用会看到价值函数在网格上出现锯齿状振荡。这也是为什么必须用迎风差分，而不是随便套个差分公式。

2.2 网格上的空间离散：左右差分和 Lax-Friedrichs 人工耗散

库内求解器处理空间导数时，会按每个坐标轴分别计算左差分和右差分。假设网格步长 dx，V 在网格点 (i, j) 上的左差分是 (V_ij - V_{i-1,j}) / dx，右差分是 (V_{i+1,j} - V_ij) / dx。y 方向同理。迎风的核心在于：梯度方向不同，采用的差分方向不同，这样信息流才符合特征传播方向。

不过只做迎风选择还不够，min/max 结构会让数值格式在非光滑点失稳，所以 ODE 求解器里普遍加一层 Lax-Friedrichs 人工耗散。下面是库里核心近似算法的简化版本：

# lf_hamiltonian.py # Lax-Friedrichs 格式的核心：给非光滑 Hamiltonian 加人工耗散 def lf_hamiltonian(H_mid, grad_plus, grad_minus, alpha): # H_mid: 用中心差分算出的 Hamiltonian 值 # grad_plus / grad_minus: 同一坐标轴上的左右差分梯度 # alpha: 人工耗散系数，取该坐标方向上状态速度的上界 return H_mid - 0.5 * alpha * (grad_plus - grad_minus)

这段代码的逻辑不复杂，但 alpha 的取值直接影响结果质量。alpha 取系统速度上界 vmax 时，耗散刚好能压住 min/max 带来的非光滑振荡；alpha 取小了，价值函数会出现波纹；取大了，价值函数会被磨成“平底锅”，捕获区域的边界模糊一大圈。我拿到一个新场景，第一件事就是把 alpha 记下来，调参时先保证它和 vmax 匹配。

2.3 时间推进与 CFL 条件：为什么步长必须跟着空间步长走

空间离散是迎风差分，时间推进用显式格式，那就绕不开 CFL 条件。库内求解器的时间步长不是固定值，而是按 dt = cfl * dx / vmax 动态算的。为什么必须这样？显式格式里，一个时间步内信息最多传播一个网格点，如果 dt 太大，信息跨过了网格，数值解就直接发散。

下面是我常用的一段参考实现，和库内逻辑一致：

# time_advance.py —— 显式时间推进与 CFL 计算 def advance(V, dt, dVdt): # V: 当前价值函数网格; dt: 时间步长; dVdt: 方程右端项 return V - dt * dVdt # 一阶显式 Euler，库内核里可换成 RK3 def cfl_dt(dx, vmax, cfl=0.5): # dx: 空间步长; vmax: 系统最大速度; cfl: 安全系数，一般取 0.3~0.5 return cfl * dx / vmax

注意 cfl 取 0.5 是保守习惯，网格加密后 dx 变小，dt 会自动跟着缩。很多人只调网格分辨率不调 dt，结果网格从 100×100 加密到 200×200 后反而炸了，原因就是 dt 没跟着 dx 一起缩。库内默认值一般安全，但你自己写扫描脚本时，每改一次网格尺寸都要确认 dt 配置是自适应的。

3. 把追逃场景跑通：编译、配置参数和读结果图的完整流程

理论看完了，接下来动手。这一章按我实际复现的路径走一遍：仓库拿到手后先看结构、再配置编译选项，然后定义一个二维追逃场景，最后把价值函数导出成图。

3.1 仓库结构与编译选项：从 C++ 内核到 Python 绑定

按我拿到的版本，常见组织方式是这样的，diffgames 用 CMake 管理，核心求解器在 src 目录，Python 绑定在 python 目录，示例场景在 examples 目录。

编译命令按标准 CMake 流程走：

# 进入 difgames 根目录后，先看 README 确认依赖版本 cmake -B build -DCMAKE_BUILD_TYPE=Release -DPYTHON_BINDINGS=ON cmake --build build -j4 # 编译完成后 Python 绑定在 build/python 下，需要把它加进 PYTHONPATH

-DCMAKE_BUILD_TYPE=Release 必须开，Debug 模式下求解迭代慢五倍以上，价值函数网格大了以后差距非常明显。-DPYTHON_BINDINGS=ON 是给后面导出结果用的，如果你只跑 C++ 示例可以关掉，但建议开着，因为后面画图、验证都依赖 Python 侧接口。编译遇到 Eigen 版本问题的话，大概率是系统装了旧版 Eigen，库内 CMake 找不到新版头文件，优先用 README 里指定的版本重新装一遍。

3.2 配置一个 200×200 的追逃场景：所有参数都在这里

我复现的第一个场景是经典的二维追逃：追逐者速度 1.0，逃逸者速度 0.8，速度比 1.25，这个配置下追逐者有理论上的捕获优势。网格范围取 [-5, 5] × [-5, 5]，节点 200×200，空间步长 dx = 0.05。

配置写在一个 YAML 文件里，大致长这样：

# scenario_two_player.yaml grid: bounds: [-5.0, 5.0, -5.0, 5.0] nodes: [200, 200] dynamics: pursuer_speed: 1.0 evader_speed: 0.8 target: radius: 0.2 numerics: cfl: 0.5 tol: 1e-6 max_iter: 5000

这里每个参数都不是随便填的。target.radius 是捕获半径，表示追逐者进入逃逸者周围 0.2 距离内就视为捕获，这个值必须大于网格步长，否则目标集在网格上可能不连通。numerics.tol 是迭代收敛阈值，看价值函数两轮迭代之间的最大变化量，小于 1e-6 就停。max_iter 设 5000 是给一个上限，正常 200×200 网格几百步就收敛了，如果一直跑满 5000 步，基本可以断定某个参数配错了。

速度比这个参数特别值得说。速度比 1.25 不是越大越好，速度比过大时捕获区域形状变化会变得很陡，网格分辨率不够会出现边界锯齿。先用 1.25 跑通，再逐步往上加，这是最稳的路径。

3.3 跑完怎么读结果：价值函数、捕获边界和收敛曲线

求解完成后导出价值函数网格：

# 求解完成后把价值函数网格导出为 npy 格式，方便用 matplotlib 画 python python/export_value_function.py results/two_player_pursuit.npy # 也可以导出为 vtk，用 ParaView 看三维曲面 python python/export_value_function.py results/two_player_pursuit.vtk

V(x) 的数值含义是捕获时间。在逃逸者周围半径 0.2 的圆内，V 等于 0，这是目标集；往外走，V 逐渐增大。画等高线图时，看 V = 1、V = 2、V = 5 这几条等高线的形状，正常情况下它们应该是围绕目标集的闭曲线；如果等高线在某个方向上开口，说明那个方向上追方没有捕获能力，对应速度比小于 1 的区域。

第一次跑完先别急着分析博弈策略，先看收敛输出里有没有报“max_iter reached”。如果迭代步数顶到上限，缩小 dt 或者检查 cfl 参数；如果求解过程出现 NaN，基本是边界条件问题，把网格边界改成外推边界再跑。输出结果和预计算的 data 目录对拍一下，V 的数值量级差在 10% 以内算正常。

4. 往场景里塞自己的规则：障碍物掩码、多智能体与速度比扫描

跑通默认场景只是开始，实际项目里要处理的是带障碍物的环境、多个智能体、还有一堆需要批量扫描的参数。这一章讲我常用的三种扩展方式。

4.1 把障碍物塞进网格：SDF 掩码和数值不可达区

带障碍物的追逃是实际项目里最常见的需求。处理方式不复杂：在价值函数网格上把障碍物区域标记为不可达。我用一个基于符号距离函数（SDF）的掩码实现：

# obstacle_mask.py —— 把障碍物写成价值函数掩码 def add_obstacle(V, grid, center, radius): # V: 当前价值函数网格; grid: 形状为 (ny, nx, 2) 的坐标网格 # center: 障碍物圆心坐标; radius: 障碍物半径 dist = np.sqrt((grid[..., 0] - center[0])**2 + (grid[..., 1] - center[1])**2) V[dist <= radius] = np.inf # 障碍物内部设为不可达，不允许进入 return V

有两点要注意。第一点，障碍物区域标记为 np.inf 而不是某个大数，因为如果是大数，梯度方向会指向障碍物内部，轨迹线反而会被吸进去；inf 让梯度在这里失去定义，轨迹自然绕开。第二点，障碍物半径至少要覆盖两到三个网格点，半径 0.05 的障碍物在 dx=0.05 的网格上只占一个点，价值函数会直接穿透它，这属于离散化精度问题，不是求解器 bug。

加完障碍物后重新求解，价值函数的等高线会把障碍物区域“挖”掉一块。验证结果是否正确有一个土办法：把最优轨迹投影到图上，看轨迹是否贴着障碍物边界走。如果轨迹离边界明显留出一大圈空隙，说明耗散系数 alpha 偏大，价值函数被磨平了；如果轨迹穿过了障碍物，说明掩码没生效，回去查 mask 的索引方向是不是反了。

4.2 多智能体不是“多一次求解”：共享价值函数与主从决策

多智能体场景最容易踩的坑是把它当成“对每个智能体各解一次方程”。实际上多智能体微分博弈里，每个智能体都要考虑对手的反应，分开求解等于把耦合丢掉了。difgames 这类开源库常见的处理方式是共享一张价值函数网格：把所有追逐者合并成一个价值函数的输入，取它们各自 Hamiltonian 的最小值，因为只要任何一个追逐者能捕获，该点就属于捕获区域。

我实际做多智能体场景时的简化做法是：如果只有一个逃逸者，多个追逐者共享同一张 V 网格，状态空间不用扩展。维度爆炸在微分博弈里是常态，双追单逃如果按完整联合状态空间算，网格维度直接翻倍，算到你怀疑人生。而共享价值函数的近似，在工程精度下损失可以接受。

另一种常见场景是主从对抗，一追一逃再加一个静态守卫者。这种我会把守卫者当作障碍物处理，不单独建博弈方程，先把主博弈算收敛，再加入守卫者的影响范围做二次掩码。它不是严格意义上的最优解，但工程上可落地，跑起来也快得多。

多智能体场景改完配置后，第一件事是画每个智能体视角下的价值函数切片。两个智能体共享一张 V 网格时，切片应该看起来相似；如果两张切片差异很大，说明你实际上在跑独立方程，耦合已经丢了。

4.3 参数扫描：速度比对捕获区域的影响

做论文配图或者方案选型时，经常要扫一组参数看趋势。我扫得最多的是速度比，因为它直接决定博弈结果的性质：速度比小于 1，捕获区域是有限闭合区域；速度比大于等于 1，捕获区域可能变成全局。批量扫描用 Python 脚本循环调用求解器：

# scan_ratio.py —— 速度比扫描 ratios = [1.05, 1.1, 1.2, 1.3, 1.5] results = {} for r in ratios: # 追逐者速度固定为 1.0，逃逸者速度取 1.0 / r V = solve(pursuer_speed=1.0, evader_speed=1.0 / r) results[r] = V

扫描完把每个速度比对应的捕获区域面积画成曲线，能看到明显的拐点：速度比 1.0 附近面积陡增，超过 1.2 后增长变缓。这类曲线放到报告里很有说服力，比贴一堆公式直观得多。

扫描最容易翻车的地方是数值参数没跟着场景变。速度比变大时，系统最大速度 vmax 也变了，如果 dt 还是按旧 vmax 算，CFL 条件可能被破坏。我的习惯是扫描脚本里每轮先重新计算 vmax 和 dt，再喂给求解器，不偷懒。

5. 避坑指南：编译、收敛和接口上最容易翻车的五个现场

这个库我用下来总体顺手，但坑也不少。以下五条都是我自己或同事实际踩过、并且能稳定复现的问题，按现象、原因、解决三段写，你照着排查省不少时间。

5.1 编译与环境问题

现象：cmake --build时直接报错，提示找不到 Eigen3 头文件，但系统里明明装了 Eigen。

原因：系统里装的是旧版 Eigen，路径和库内 CMake 要找的版本对不上。Eigen 是头文件库，版本新旧不体现在.so 文件上，只体现在头文件路径和宏定义里，CMake 找到旧路径后不会主动报版本不兼容，等到编译某个用新 API 的源文件时才炸。

解决：按 README 指定的 Eigen 版本重新安装，并在 cmake 命令里显式指定路径，-DEIGEN3_INCLUDE_DIR=/path/to/eigen3，不要再依赖系统自动查找。装完先编译一个空示例验证头文件路径生效，再编主项目。

5.2 数值与收敛问题

现象：网格从 100×100 加密到 200×200 后，价值函数反而出现明显的锯齿振荡，迭代步数顶到 max_iter 也不收敛。

原因：网格加密后 dx 变小，但 dt 还是按旧网格算的。显式时间推进的 CFL 条件被破坏，误差在每个时间步累积，最终在非光滑 Hamiltonian 附近激发振荡。

解决：网格尺寸一改，dx、dt、alpha 三个参数必须同步更新。我一般把 dt 和 alpha 都做成网格参数的函数，写在配置文件里而不是硬编码，避免手改网格时漏掉。

现象：目标集周围的价值函数出现负值，看起来像“凹下去”的不自然形状。

原因：初始值处理不对。价值函数初始化成全零时，目标集外区域的梯度信息缺失，迭代过程会把目标集附近的 V 推成负值。库内求解器对初值有默认处理，但我自己写扩展场景时经常踩。

解决：把初始值设置成到目标集的符号距离函数，而不是全零。这样迭代一开始梯度就有正确的量级，目标集附近不会产生非物理负值。

5.3 场景与接口问题

现象：半径 0.03 的小障碍物完全没起作用，轨迹直接穿过。

原因：dx = 0.05 时，半径 0.03 的圆形障碍物在网格上覆盖不了任何网格点，掩码数组全为 False，等于没加障碍物。这不是求解器 bug，是离散化精度问题。

解决：障碍物半径至少取网格步长的 2 到 3 倍；小于这个值要么缩小 dx，要么把障碍物当作软约束按大数惩罚处理，不要硬塞 mask。

现象：Python 绑定里开了多个线程同时调求解器，程序直接卡死，CPU 占用率上不去。

原因：Python 绑定内部持有 GIL，多线程调用求解器时实际上串行执行，并且线程切换频繁导致死锁概率大增。

解决：多场景并行改用 multiprocessing，每个进程独立持有解释器；或者把批量求解循环全部放进 C++ 侧，Python 只负责发起和收集结果。

6. 验证价值函数正确性：一条不依赖仿真器的轨迹回溯法

价值函数算完，怎么确认它是对的？你当然可以写一个完整的博弈仿真器来验证，但那样成本太高。我常用的是一个轻量级回溯法：对着收敛后的 V 网格，从任意起点沿负梯度方向做数值积分，得到一条最优轨迹；然后交叉验证两点——轨迹是否落在目标集上，以及轨迹累计时间是否约等于 V(x0)。

# verify_value_function.py —— 轨迹回溯验证 def verify_trajectory(V, x0, target_set, params): # V: 已收敛的价值函数网格; x0: 起始位置 # target_set: 目标集掩码; params: 包含网格步长和时间步长 x = np.array(x0, dtype=float) path = [x.copy()] for _ in range(int(2.0 / params["dt"])): # 判断当前位置是否进入目标集 iy = int((x[1] - params["bounds"][0]) / params["dx"]) ix = int((x[0] - params["bounds"][1]) / params["dx"]) if target_set[iy, ix]: break # 用网格梯度反推最优控制方向，取负梯度并归一化 gy, gx = np.gradient(V, params["dx"], params["dy"]) g = np.array([gx[iy, ix], gy[iy, ix]]) if np.linalg.norm(g) < 1e-9: break # 梯度为零，说明该点可能不在价值函数有效域内 direction = -g / np.linalg.norm(g) x += params["dt"] * params["vmax"] * direction path.append(x.copy()) return np.array(path)

逻辑很直接：V 的梯度方向就是最优控制方向，从任一点出发不断沿负梯度走，最后应该落到目标集边界。如果轨迹终点离目标集很远，说明价值函数在某个区域的梯度指向错了，优先回查该区域的网格分辨率和耗散系数；如果轨迹长度乘以时间步长和 V(x0) 差超过百分之十，大概率是耗散系数偏大，价值函数被磨平，导致路径时间短于理论捕获时间。这个方法不需要博弈仿真器，一套对比结果就能定位问题出在离散还是出在参数。

从那以后我每次换场景、调参数，都会先把这套回溯验证跑一遍，再谈别的。V(x0) 和轨迹时间对不上，值再好看我也不会拿去用。希望帮到你。

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