简介:微分博弈开源项目聚焦动态系统中的多人策略决策问题,以哈密顿-雅可比-贝尔曼-伊萨克斯(HJBI)方程为理论核心,提供适用于Linux平台的C语言实现。代码中包含多种数值求解与模拟场景,用于求解连续时间博弈的纳什均衡;配套图形模块可直观展示博弈进程,Shell脚本则支持一键编译运行,PDF文档说明了优化原理与实验方法,其中数值求解可能涉及有限差分、有限元等近似手段。整个压缩包共9个文件,以C源文件为主(6个.c),并含Shell脚本、PDF文档和头文件,总大小约206KB。项目结构紧凑,适合博弈论、控制理论与计算数学方向的开发者或研究者阅读源码、复现实验,并可将相关解法迁移至自动驾驶、经济建模等实际场景。目前已有238人学习浏览,从代码到文档一应俱全,是理解HJBI方程工程实现的实用参考资料。
1. 微分游戏的开源实现:difgames 把 HJI 方程变成能跑的追逃博弈仿真
微分游戏(differential game)听起来像数学课上的抽象名词,但第一次在仿真里看到追逃双方的轨迹贴着目标集边界走,我才意识到这套理论其实非常工程化。difgames 是一个聚焦微分博弈的开源软件框架,核心只做两件事:在网格上数值求解 HJI 方程,再把值函数转成可达集、安全集与最优控制策略。它适合正在做机器人避障、自动驾驶决策、多智能体对抗仿真的人——你不必从零写偏微分方程求解器,也不用把论文里的公式一行行手工翻译成代码,就能复现出追逃博弈里的控制律。下文按“数学内核 → 快速上手 → 参数调优 → 避坑 → 验证”推进,跟着走可以跑通第一个二维追逃实例。
2. 数学内核与代码结构:HJI 方程在网格上是怎样被演化出来的
2.1 零和博弈与 HJI 方程:difgames 到底在解什么
追逃博弈是标准的零和博弈:追捕者想让终端时刻两方距离尽可能小,逃逸者想让距离尽可能大。这个“一个想最小、一个想最大”的对抗关系,在数学上被写成值函数 V(x,t),它表示从状态 x、时刻 t 出发,在双方都采取最优策略后的终端代价。difgames 把求解这个值函数作为一切功能的地基,得到 V 之后,可达集、安全集、最优控制输入都是对 V 的后续处理。
HJI 方程是这个框架真正在解的偏微分方程,形式上是时间反向的哈密顿-雅可比方程。对零和博弈,方程里会出现一个 min 和一个 max,分别对应两个玩家的优化方向。工程上最容易搞混的是 min 和 max 的顺序:顺序写反,就是让逃逸者去最小化、追捕者去最大化,解出来的值函数会变成另一种博弈的结果,而且从数值曲面上看很难一眼识别。
difgames 里用solver_type参数区分零和与非零和博弈,默认走零和分支。实际项目中,我一般会先在注释里把两个玩家的角色写清楚,再对照源码确认u和d两个输入分别被放在 min 还是 max 一侧。因为这类框架的接口通常以“控制”和“扰动”来命名,追捕者对应控制,逃逸者对应扰动,一旦把接线接反,后面所有分析都建立在错误的值函数之上。
方程只在网格点上被离散求解,这意味着输入的状态空间必须是有界的。difgames 要求每个维度都给定lo和hi,超出边界的状态由边界条件处理。物理上这对应“状态范围假设”,比如位置在正负两米内、速度在正负一米每秒内。边界设小了会把真实可达集切掉,设大了会在不需要的地方浪费算力,这一步往往是新手第一个翻车点。
2.2 水平集方法与 LF 格式:值函数演化是一张曲面在时间轴上倒放
目标集在 difgames 里被表示成一个函数 g(x),比如“到原点距离小于 0.3”就是 g(x)=||x||−0.3。g(x)≤0 的点构成目标集,g(x)=0 的等值面就是目标集边界。这个函数同时还是值函数在终端时刻的初始条件,求解过程就是让这张曲面在时间轴上从终端时刻向当前时刻反向演化。
这种用等值面隐式表达几何的方式,就是常说的水平集方法。它的好处是不用显式追踪边界上的点,边界怎么拓扑变化都不需要额外处理,天然适合处理两个可达集合并、撕裂这类情况。代价是必须有一个足够密的网格,曲面细节小于网格间距时会被磨平,导致目标集“消失”,这一点后面避坑章节会专门讲。
时间方向上的演化用的是显式格式,空间梯度用中心差分近似,再叠加一个人工扩散项保证数值稳定。这个组合在行业里常被叫做 LF 格式,是求解 HJ 方程最省事的方案,虽然精度只有一阶,但胜在实现简单、对内存访问友好,而且天然适合 Numba 或 C++ 加速。difgames 默认后端就是 Numba,纯 NumPy 后端用于调试。
从工程视角看,LF 格式的人工扩散系数给得越大越稳定,但解会被抹得更平,边界不再锐利。我一般会把这个系数取成与状态导数上界同一量级,而不是按默认值一跑了之。经验是:如果你看到等值面边缘出现一条条细密的平行波纹,通常是人工扩散太小或空间步长太大,而不是模型本身的问题。
2.3 三个核心类:Grid、Dynamics 与 Solver 的分工
打开 difgames 源码,核心抽象就三个类:Grid描述网格,Dynamics描述状态方程,Solver负责把两者组合起来做时间演化。理解这三个类的分工,基本就理解了这个框架的扩展方式。
# core_interfaces.py(difgames 的类接口风格示意) from abc import ABC, abstractmethod import numpy as np class Dynamics(ABC): def __init__(self, state_dim, control_dim, disturbance_dim): self.state_dim = state_dim self.control_dim = control_dim self.disturbance_dim = disturbance_dim @abstractmethod def dxdt(self, x, u, d): """返回状态导数,u 是己方控制,d 是对手/扰动输入。""" class Grid: def __init__(self, lo, hi, n): self.lo = np.asarray(lo, dtype=float) self.hi = np.asarray(hi, dtype=float) self.n = np.asarray(n, dtype=int) self.dx = (self.hi - self.lo) / self.n self.axes = [np.linspace(self.lo[i], self.hi[i], self.n[i]) for i in range(len(self.lo))] def shape(self): return tuple(self.n) class Solver: def __init__(self, grid, dyn, t_max, dt, backend="numba"): self.grid = grid self.dyn = dyn self.t_max = t_max self.dt = dt self.backend = backend def solve(self, value_init): """从终端时刻的值函数出发,反向演化到当前时刻。""" value = np.asarray(value_init) steps = int(round(self.t_max / self.dt)) for _ in range(steps): value = self._step(value) return value代码里的value_init就是上文的 g(x),通常在业务代码里用一个 lambda 或独立函数构造出来。_step在真实源码里会调用 LF 格式的核函数,并对边界做填充;我这里只留了接口占位,实际项目里你会看到它内部还有_hamiltonian和_boundary_extend两个私有方法,分别负责组装带符号的数值通量和处理边界外延。
这里要注意:solve返回的是value,它的语义是“从当前状态出发、在这个时间窗口内能到达目标集的可能性”。严格说是值越小越容易到达,值为负表示该状态已经在可达集内部。很多第一次用的人会把值函数和代价函数搞混,以为值越大越好,这是理解后续所有后处理的一个隐性门槛。
3. 快速上手:搭一个二维追逃博弈仿真实战
3.1 环境准备与安装
difgames 的源码包到手后,第一步是建虚拟环境再以可编辑模式安装。可编辑模式对这类数值框架特别重要,因为你大概率会改到 Dynamics 子类或数值核函数,每次改完不用重新 install,Python 直接就能载入最新代码。
# 进入 difgames 源码根目录后执行 python -m venv .venv source .venv/bin/activate # Windows 下执行 .venv\Scripts\activate pip install --upgrade pip pip install -e .安装时 pip 会自动拉起 numpy、scipy、numba 三个核心依赖。numba 是默认后端,如果机器上装 LLVM 有困难,可以在环境变量里指定退回到 NumPy 后端:export DIFGAMES_BACKEND=numpy。NumPy 后端慢大约一个量级,但跑 2 维小网格调试是完全够用的。
3.2 定义双积分器动力学模型
二维双积分器是最常见的调试模型:状态是位置和速度四个维度,控制是 x 和 y 方向的加速度。用它验证算法,既能看到轨迹的惯性特征,又不至于像完整车辆模型那样有一堆轮胎力参数要调。
# dynamics/double_integrator.py import numpy as np from difgames.core import Dynamics class DoubleIntegrator(Dynamics): def __init__(self, acc_bound=1.0): super().__init__(state_dim=4, control_dim=2, disturbance_dim=2) self.acc_bound = acc_bound def dxdt(self, x, u, d=None): px, py, vx, vy = x ax, ay = np.clip(u, -self.acc_bound, self.acc_bound) if d is not None: ax += d[0] ay += d[1] return np.array([vx, vy, ax, ay])state_dim=4是位置加速度四个自由度,disturbance_dim=2给逃逸者留了两个加速度通道。真实追逃场景里追捕者控制加速度,逃逸者的策略则通过d注入到同一个动力学上,这样 HJI 方程里 min 和 max 作用在同一组状态方程上,物理意义才是同一个追逃系统,而不是两套独立轨迹。
np.clip的边界直接写在dxdt里,是为了保证控制量在任何后端下都合法。数值求解 HJI 方程时,Hamiltonian 里也会用到控制界,两处必须保持一致,否则会出现可控集计算和轨迹仿真对不上的问题。这个参数在acc_bound上只设一处,后续所有代码都从self.acc_bound读取,避免两处手改漏掉一处。
3.3 配置目标集与启动仿真
追逃场景里目标集通常是“逃逸者初始所在位置的邻域”。为了让结果直观,我常把目标集设在原点附近的一个小球,追捕者从远处出发,看看在指定时间窗内能覆盖多大范围。
# run_pe.py import numpy as np from difgames.core import Grid, Solver from dynamics.double_integrator import DoubleIntegrator grid = Grid(lo=[-2.0, -2.0, -1.0, -1.0], hi=[ 2.0, 2.0, 1.0, 1.0], n=[80, 80, 60, 60]) dyn = DoubleIntegrator(acc_bound=1.0) def target_set(x): px, py, vx, vy = x return np.sqrt(px**2 + py**2) - 0.3 solver = Solver(grid, dyn, t_max=1.5, dt=0.01) value = solver.solve(target_set) # 提取可达集边界:值是 0 的等值面 print("min value:", value.min(), "max value:", value.max()) print("reachable set cell count:", int((value <= 0).sum()))target_set用欧几里得距离减半径的方式定义,返回值小于等于零表示在目标集内部。四个维度的网格中,位置维度范围正负 2 米,速度维度范围正负 1 米每秒,速度维度点数少一些,因为速度方向上值函数的变化通常比位置方向平缓。
value <= 0的网格点数量就是可达体积的粗略估计。这一步先不看具体数值,只看它是否在一个合理量级——如果算出来是 0 或者全网格点都小于 0,说明参数配错了,应该回到目标集和网格配置上排查。常见问题是你以为在跑追捕者视角,实际上逃逸者视角下的可达集形状完全不同。
4. 参数怎么设:网格分辨率、时间步长与终端时间
4.1 网格分辨率:算力与精度的第一道坎
网格点数直接决定内存和耗时,公式很简单:存储一个 float64 值函数需要 8 字节乘以网格总点数。但实际计算时并不只存一份数组,LF 格式至少要同时保留当前层、下一层和梯度计算用的临时数组,所以内存占用通常是单份数组的四到六倍。
按前面那个 80×80×60×60 的网格算,总点数是 2304 万,单份数组约 184 MB,乘上四倍就是超过 700 MB。这已经是一个很吃内存的配置了。我一般调试时先从 40×40×30×30 起步,确认值和形状合理后,再翻倍到 60×60×40×40,最后才跑完整分辨率。
还有一个反直觉的点:速度维度的分辨率对追逃可达集的影响比位置维度更大。因为值函数在速度方向上的梯度往往更大,网格太粗会把速度边界磨掉,导致可达集在速度方向上“缺一块”。如果内存有限,优先保速度维度点数,位置维度可以适当缩。
4.2 时间步长:CFL 条件与数值稳定性
显式时间格式都有稳定性约束,LF 格式的约束表现为时间步长和空间步长必须满足类似 CFL 的条件。粗略估计方法是用最小网格间距除以状态导数的最大模,得到一个步长上界,再乘上 0.8 的安全系数。这个“状态导数最大模”在双积分器里大致是速度上限加加速度上限乘以时间步长,但加速度项本身就是隐式耦合,所以我会再留 20% 余量。
实操中我常用一个更保守的公式:dt = 0.5 * min(grid.dx) / v_max,其中v_max是状态空间里能出现的最大速度。对双积分器,v_max就是速度维度的上界;对带扰动的系统,还要看扰动对速度影响的最大幅值。跑出来的结果如果出现 NaN,第一步不是怀疑模型,而是把 dt 再减半试一次。
注意 dt 也不是越小越好,时间步长小了,计算步数线性增加,而 LF 格式的人工扩散误差会随着迭代步数累积。最佳实践是先满足稳定性下界,然后同时增大空间分辨率,而不是单独缩小 dt。因为只有同时缩小 dx 和 dt,才能保证真实解的误差在下降。
4.3 终端时间与收敛判据
终端时间 t_max 决定了可达集能有多大。双积分器从静止出发一秒内能移动 0.5 米,但如果有初速度,同样一秒能移动更远。理论上时间越长可达集越大,但增长会趋于饱和,因为控制输入有限,系统最终会被速度边界卡住,值函数不再随时间变化。
我判断 t_max 是否足够的办法是跑两组对比:分别用 t_max 和 2 倍 t_max 求解,观察整个网格上值函数的最大变化量。如果相差小于 1e-6,说明已经收敛到稳态解,再增加时间也不会让可达集明显扩大。如果差距很大,说明时间窗口还处于瞬态阶段,输出结果对终端时间很敏感,这时候直接拿去画可达集是危险的。
下表整理了我在 difgames 里常用的一组起点参数,适合四维双积分器模型,后面再根据自己的算力逐步放大。
| 参数 | 建议起点 | 说明 |
|---|---|---|
| 位置网格 n | 60 每维 | 位置方向等势面较平滑 |
| 速度网格 n | 40 每维 | 速度方向梯度大,优先保点数 |
| dt | 0.5 * min(dx) / v_max | 乘 0.8 安全系数 |
| t_max | 收敛判断后确定 | 对比 t_max 与 2*t_max 的差 |
| 目标集半径 | 至少 3 * dx | 太小会被网格采样漏掉 |
5. 避坑与排查:五条能省下一周排错的实战经验
5.1 值函数发散成 NaN,残差越跑越大
现象:solve()跑完,value数组里出现大片 NaN,或者整个值变成正负无穷。
原因:时间步长超过了 LF 格式的稳定域,显式时间推进产生数值振荡,振幅逐层放大直到溢出。我在双积分器模型上第一次遇到时,还怀疑是动力学函数里有除零,后来检查dxdt发现干净得很,其实就是 dt 设得太大。
解决:把 dt 降到0.5 * min(dx) / v_max的 60% 再跑。如果还不稳定,检查target_set是否有奇异值,目标函数在网格边界上出现无穷大也会把整层值函数带崩。可以在solve循环里加一个np.isfinite(value).all()断言,直接在第一次溢出时告诉你第几步炸的。
5.2 目标集边界出现锯齿状波纹,等值面不光滑
现象:画出来的可达集边界不是一条平滑曲线,而是像毛边一样参差不齐,但值函数本身没有 NaN。
原因:人工扩散系数偏小,中心差分格式没有足够阻尼,高频分量没有被压下去。另一个常见原因是网格各方向间距差异过大,比如位置方向 dx=0.05 而速度方向 dx=0.02,LF 格式的各向异性扩散会优先抹平速度方向,导致边界在位置方向上留下毛刺。
解决:把人工扩散系数从默认值调大 30% 到 50%,观察波纹是否消失。注意不要为了磨平边界把系数加太多,否则整个目标集会膨胀一圈,看起来光滑但已经是错误解。我一般会做一组系数扫描,选“刚好消除波纹”的最小值,这是最不损害精度的方案。
5.3 min 和 max 顺序写反,追捕者和逃逸者角色颠倒
现象:仿真轨迹里逃逸者主动冲向追捕者,或者可达集的形状和预期完全反向——该大的地方小,该小的地方大。
原因:Dynamics.dxdt里u和d的符号被接反,HJI 方程的 Hamilton 项里原本应该取 min 的变成了取 max。这种错误在数值上不会报错,甚至在二维图上也能画出漂亮的等值面,但物理意义完全错了,最难排查。
解决:先在最简单的一维单积分器上做单元测试:固定逃逸者策略为零,只优化追捕者控制,看值函数是否等于解析可达集。跑通后再加入对抗项。我现在的习惯是在run_pe.py顶部用注释明确写出“u 是追捕者加速度,d 是逃逸者加速度”,防止三个月后自己回头看也忘了当初的设定。
5.4 目标集半径太小,整个值函数恒为零
现象:value.min()和value.max()都接近 0,可达集计数为 0,好像系统根本没有任何控制能力。
原因:目标集半径小于网格间距,初始化时没有任何网格点落在目标集内部,值函数在所有网格点上都是正的。水平集方法用网格点上的采样值代表连续曲面,目标集比网格粒度还小的时候,它就“看不见”这个目标了。
解决:把目标集半径调大到至少三倍网格间距,比如 dx=0.05 时半径取 0.15。如果需要精确表达更小的目标,就不能靠加密网格硬扛,而是改用连续符号距离函数做初始化,让等值面穿过网格点之间而不是恰好落在空白处。
5.5 内存不够,一跑就 OOM
现象:网格设置稍大,比如 100×100×60×60,进程直接内存溢出被杀掉,连报错都来不及。
原因:前面提到的多份临时数组,加上 NumPy 每个网格节点会分配对象头部开销,实际内存占用常常是预期的五倍以上。我在一台内存不大的机器上跑过 100×100×60×60,单份数组 288 MB,进程峰值直接冲到 1.6 GB,然后被杀。
解决:先把后端切到 NumPy 调试模式,用小网格跑通逻辑,再切 Numba 后端跑大网格。同时把不用的临时数组用del及时释放,或者把网格分成块来求解。Numba 后端还有自己的内存复用机制,开启后峰值内存能降一半,但代价是 JIT 编译时需要额外的编译时间。
6. 验证方法:用单积分器解析解给数值结果做体检
在把 difgames 用到复杂模型之前,我会强制跑一遍单积分器的解析验证。单积分器状态只有 x、y 两个维度,控制是速度且最大值 1,目标集是半径 r0=0.2 的圆。解析上从任意点出发在 T 秒内能到达目标集的充要条件是初始位置到原点距离不超过 r0+T,因此终端时刻的值函数应当是 V(x)=||x||−(r0+T)。
import numpy as np from difgames.core import Grid, Solver class SingleIntegrator(Dynamics): def dxdt(self, x, u, d=None): return u def verify(T=1.0, r0=0.2, n=100): grid = Grid(lo=[-2, -2], hi=[2, 2], n=[n, n]) dyn = SingleIntegrator(state_dim=2, control_dim=2, disturbance_dim=0) vmax = 1.0 solver = Solver(grid, dyn, t_max=T, dt=0.5 * grid.dx[0] / vmax) value = solver.solve(lambda x: np.sqrt(x[0]**2 + x[1]**2) - r0) mesh0, mesh1 = np.meshgrid(grid.axes[0], grid.axes[1], indexing="ij") analytic = np.sqrt(mesh0**2 + mesh1**2) - (r0 + T) return np.mean(np.abs(value - analytic)) for n in [40, 80, 160, 320]: err = verify(n=n) print(f"n={n}, mean_abs_err={err:.4f}")这个测试跑出的误差应当随网格加密而下降,而且每翻一倍网格点数,误差大致减半,对应 LF 格式一阶收敛的特征。如果误差不降甚至升高,说明边界条件、人工扩散系数或符号方向有一个是有问题的。我第一次跑这个验证时,误差在 160 格后不降反升,检查半天发现是边界填充方式用了固定值,把真实可达集反射回了网格内部。
从那以后我每次拿到新版本 difgames,或者换一套动力学模型,都要先跑一遍这个圆形验证再动业务参数。否则在复杂模型里遇到反直觉的结果,你根本分不清是参数没调好还是框架用法有误,而这个最小复现用例能帮你快速定位问题到底出在哪个模块。希望帮到你。
本文还有配套的精品资源,点击获取