简介这是一份面向博弈论、控制理论及计算数学研究者的开源微分博弈项目。项目以哈密顿-雅可比-贝尔曼-伊萨克斯HJBI方程为核心展示如何用数值方法求解动态博弈中的最优策略与纳什均衡适合研究生、算法工程师及对自动驾驶、经济建模等应用场景感兴趣的开发者学习参考。资源共含9个文件以C语言源码为主覆盖求解算法实现、图形可视化模块与文献管理工具另有PDF文档说明优化方法与实验结果并提供build.sh脚本帮助用户在Linux环境下快速编译运行。压缩包仅206KB代码结构紧凑便于直接阅读与二次开发。目前已有237人学习下载可作为一个轻量而完整的入门样例帮助读者理解HJBI方程的建模思路、数值求解流程及动态博弈结果的图形化展示方法。1. 微分博弈的求解门槛difgames 这个开源库替你做了哪件事做机器人对抗、无人机追逃方向的工程师早晚会撞上一堵墙博弈逻辑好写但 Hamilton-Jacobi-IsaacsHJI方程数值解不出来。这个方程里带 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 人工耗散库内求解器处理空间导数时会按每个坐标轴分别计算左差分和右差分。假设网格步长 dxV 在网格点 (i, j) 上的左差分是 (V_ij - V_{i-1,j}) / dx右差分是 (V_{i1,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, cfl0.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 目录。目录内容我一般怎么用src/C 数值内核网格与 HJI 求解器不直接改看清楚接口就行examples/二维追逃、避障、多智能体示例改参数最常从这里开始python/Python 绑定与导出脚本画图、跑批处理都走这里data/预计算的价值函数结果先看结果再跑自己场景能对拍编译命令按标准 CMake 流程走# 进入 difgames 根目录后先看 README 确认依赖版本 cmake -B build -DCMAKE_BUILD_TYPERelease -DPYTHON_BINDINGSON cmake --build build -j4 # 编译完成后 Python 绑定在 build/python 下需要把它加进 PYTHONPATH-DCMAKE_BUILD_TYPERelease 必须开Debug 模式下求解迭代慢五倍以上价值函数网格大了以后差距非常明显。-DPYTHON_BINDINGSON 是给后面导出结果用的如果你只跑 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.vtkV(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 的障碍物在 dx0.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_speed1.0, evader_speed1.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, dtypefloat) 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) 和轨迹时间对不上值再好看我也不会拿去用。希望帮到你。本文还有配套的精品资源点击获取
