简介微分游戏difgames开源项目面向博弈论、控制理论与动态优化方向的学习者和研究人员聚焦连续时间动态系统与博弈论结合下的HJBI方程求解与仿真实践适合具备一定C语言和数学基础、希望动手理解纳什均衡与最优策略的读者。压缩包共9个文件约206KB以6个C源码为主体配合1个头文件、1个Shell构建脚本和1份PDF文档源码承载算法与图形显示脚本负责自动化编译文档则补充优化方法与实验说明。已有239人学习说明该主题在相关社区具有一定关注度。读者可从中获得HJBI方程的C语言实现思路、偏微分方程的数值近似方法、动态博弈过程的可视化呈现以及构建运行的环境入口并借此把抽象理论落到可编译、可观察的代码实验中为自动驾驶、经济模型或策略分析等场景提供参考。1. 追逃对抗算不出反馈策略先看 difgames 这套开源软件把 HJI 拆成了什么追逃这类问题单智能体最优控制那套办法不太好使对手也在优化你算出的最短路径随时会被对方的策略改写最后得到的是一个鞍点而不是一条轨迹。difgames 是一款把微分游戏differential game从公式变成可运行代码的开源软件核心逻辑只有三件事在状态空间上离散并求解 Hamilton-Jacobi-Isaacs 方程得到价值函数再从价值函数的梯度中提取双方的反馈策略。它适合两类人做机器人或无人系统对抗博弈、需要回答“到底能不能追上、多久追上”的工程同学以及课题里需要绘制可达集、需要拿解析解做校验的研究同学。下面按我拆包的顺序写HJI 如何落到网格上第一个追逃算例怎么跑通参数怎么标定以及那几个能把人熬到后半夜的坑。2. 把微分博弈落成可算的网格HJI 方程、Isaacs 条件与离散化选型2.1 谁在 min、谁在 maxHJI 的符号约定两方零和追捕方 u 想让代价最小逃逸方 v 想让代价最大动力学写成 ẋ f(x, u, v)u ∈ U、v ∈ V 是各自的可行控制集。价值函数 V(x, t) 表示从状态 x、时刻 t 出发的最优代价满足∂V/∂t min_{u∈U} max_{v∈V} { ∇V · f(x, u, v) } 0 V(x, T) g(x)这里最容易翻车的是符号和顺序。min 在外、max 在内对应“追捕方先动、逃逸方后动”反过来写成 max min只有在 Isaacs 条件成立时两者才相等也就是哈密顿量存在鞍点。工程上判断起来不难如果 f 能拆成 f₁(x, u) f₂(x, v) 这种可分离形式鞍点自动存在直接取 min max 就行如果是强耦合的比如推力方向同时受双方姿态影响就得老老实实检查鞍点条件否则解出来的根本不是博弈解。我一般的做法是把 V 定义成“追捕方还差多少代价才能完成任务”这样 V 小的地方对追捕方有利。追捕方的反馈策略沿 −∇V 走逃逸方沿 ∇V 走方向反了会让闭环仿真里的追捕器绕着目标转圈——这个现象后面第 5 章还会提到。2.2 为什么走可达集而不是直接解策略新手最容易想到的思路是既然要策略那就直接对策略做优化、做参数化用强化学习或者打靶法去搜。这条路在单项任务上能出结果但拿不到“保证”你搜出来的策略只在采样到的初始条件附近可信换一个初始位置就得重训。网格型 HJI 求解器给的是全局的。一次求解整个状态空间上每个点的价值函数都有值零水平集就是“追捕方能在时限内完成任务”的状态集合也就是后向可达集。它天然处理非光滑追逃博弈的价值函数在“屏障面”barrier上会出现折痕最优策略在那里切换用特征线法硬推会在折痕处直接崩掉而网格法 数值粘性恰好收敛到粘性解这是这套方法最值钱的地方。代价是维数灾难。二维状态网格 n201 是四万个点跑起来很轻一旦你把追捕方和逃逸方各自的二维位置都塞进网格变成四维n81 就是四千三百万个点双精度直接 340MB 起再往上就等着 swap。所以选型上的分水岭很清楚动力学是简单运动各向同性、速度恒定时把状态约化成相对位置二维就够精度和速度都能接受一旦有转弯半径约束、加速度上限这类各向异性就只能上全状态或者增广状态这时候要么降网格、要么上 GPU。2.3 环境与最小可运行骨架依赖很薄numpy 做数组、scipy 做插值和稀疏运算、matplotlib 看零水平集。源码包拿到手后常见做法是先建虚拟环境再以可编辑模式装上这样改求解器不用反复重装。python -m venv .venv source .venv/bin/activate pip install numpy scipy matplotlib # 进到解压后的项目根目录以可编辑模式安装方便直接改求解器源码 pip install -e . # 先跑最小示例确认依赖齐了、示例数据在 python examples/pursuit_simple.py --grid 101 --horizon 1.0三个参数分别对应--grid是每维网格点数先用 101 探路确认能出图再上 201--horizon是逆向推进的时长物理含义是“我关心多久之内能不能追上”示例脚本名以你手上的目录为准认准examples/下带 pursuit 或 reachability 的那个就行。跑通之后先别急着换自己的动力学。包里自带的求解器核心就是把“建网格—算哈密顿量—时间推进”这三段拼起来我把这三段单独摘出来写了一份不依赖包的等价实现用来跟库里的结果做对照后面第 3 章直接给代码。这一步看着多余实际是整个流程里最省时间的一步库的输出对不上时你至少知道自己手算的那份是对的。3. 跑通第一个追逃算例网格、终端条件与时间推进拿最简单的追逃做验证追捕方速度 a、逃逸方速度 b都取全向相对位置 r 满足 ṙ b·v − a·u|u| |v| 1。这个算例有解析解是天然的单元测试。3.1 建网格与终端条件终端条件 g(x) 用“到目标集的有符号距离”目标集内部为负外部为正。追逃问题里目标集就是原点所以 g(x) |x|零水平集退化成一个点。import numpy as np def build_grid(n, L): n: 每维点数; L: 计算域半边长 x np.linspace(-L, L, n) dx x[1] - x[0] X, Y np.meshgrid(x, x, indexingij) return X, Y, dx def terminal_value(X, Y): 终端条件: 到捕获点(原点)的距离, 零水平集 原点 return np.hypot(X, Y)有符号距离用|x|而不是|x|²很关键。平方距离在原点附近梯度趋零数值格式会在这里丢掉一阶精度零水平集被抹成一团取距离本身梯度模长恒为 1格式在原点附近反而最准。计算域半边长 L 的取法初始相对距离的 1.5 到 2 倍。太小零水平集在时限内会撞到边界边界条件会污染结果太大同样网格数下 dx 变粗界面精度又不够。3.2 哈密顿量的实现对上面的相对动力学min over u 和 max over v 可以解析算出来min_u ∇V·(−a u) −a|∇V|max_v ∇V·(b v) b|∇V|合起来H(∇V) (b − a) · |∇V|a b 时 H 为负逆向推进中零水平集以 (a−b) 的速度向外扩张——这就是“追捕方更快捕获区随时间长大”。a b 时符号翻转集合收缩对应逃逸方占优、初始不在捕获区就永远追不上。def hji_rhs(V, dx, a, b): 返回 V 对逆向时间 tau 的导数: V_tau -H 数值粘性 Vx np.zeros_like(V); Vy np.zeros_like(V) Vx[1:-1, :] (V[2:, :] - V[:-2, :]) / (2 * dx) Vy[:, 1:-1] (V[:, 2:] - V[:, :-2]) / (2 * dx) p np.hypot(Vx, Vy) # |∇V| H (b - a) * p # 可分离动力学的鞍点哈密顿量 # Lax-Friedrichs 型数值粘性, alpha 取哈密顿量对梯度的最大特征速度 alpha abs(a - b) lap np.zeros_like(V) lap[1:-1, 1:-1] (V[2:, 1:-1] V[:-2, 1:-1] V[1:-1, 2:] V[1:-1, :-2] - 4 * V[1:-1, 1:-1]) / dx**2 return -H 0.5 * alpha * dx * lap粘性项系数里的0.5 * alpha * dx不是拍脑袋定的Lax-Friedrichs 通量拆开之后耗散项的物理形式就是 (α/2)·dx·ΔVα 取哈密顿量对梯度偏导的最大模。这一步写错界面要么抖成马蜂窝要么糊成一片。3.3 时间推进与 CFL逆向时间 τ 从 0 走到 Tτ0 时是终端条件τ 越大捕获集越大。步长由 CFL 定def solve(n201, L2.0, a1.0, b0.6, T1.5, cfl0.4): X, Y, dx build_grid(n, L) V terminal_value(X, Y) speed max(abs(a - b), 1e-12) dt cfl * dx / speed steps int(np.ceil(T / dt)) dt T / steps # 重新均分, 保证正好走到 T for k in range(steps): V V dt * hji_rhs(V, dx, a, b) # 边界线性外推: 外侧比内侧多走一个 dx 的距离, 近似无界域 V[0, :] V[1, :] dx; V[-1, :] V[-2, :] dx V[:, 0] V[:, 1] dx; V[:, -1] V[:, -2] dx return X, Y, V, T每一行都值得停一下。dt T / steps这步重算不是洁癖如果直接用 CFL 给的 dt 跑循环、最后一帧停在 T 之前你拿去跟解析解对拍时误差里会混进“时间没走完”的系统偏差查半天查不出原因。边界外推成内邻 dx而不是固定成某个大常数是因为这里的价值函数在无界域上本来就以单位斜率增长外推等价于把远处的距离继续算下去。3.4 从价值函数里抽出反馈策略有了 V策略就是梯度方向def feedback(V, dx, a, b, eps1e-9): Vx np.zeros_like(V); Vy np.zeros_like(V) Vx[1:-1, :] (V[2:, :] - V[:-2, :]) / (2 * dx) Vy[:, 1:-1] (V[:, 2:] - V[:, :-2]) / (2 * dx) g np.maximum(np.hypot(Vx, Vy), eps) # 防止除零, 原点附近梯度退化 ux, uy -Vx / g, -Vy / g # 追捕方: 沿 -∇V, 让 V 变小 vx, vy Vx / g, Vy / g # 逃逸方: 沿 ∇V, 让 V 变大 return ux, uy, vx, vyeps这个保护项是必需的。零水平集上 |∇V| 理论上不为零距离函数的梯度恒为 1但数值解在原点附近会因为网格对称性出现极小值直接除会得到方向乱跳的单位向量闭环仿真里表现为追捕器在目标附近抽搐。加了下限之后靠近原点的方向精度下降但至少不会给出反向指令。跑完这三段你手上应该有一张零水平集图一个以 (a−b)τ 为半径向外扩张的圆。圆不圆、扩张速度对不对就是下一章参数标定的起点。4. 参数标定dt、CFL 系数、数值粘性与收敛判据怎么定4.1 网格分辨率与内存的取舍先说清代价。二维网格下n 翻倍内存是四倍、计算量也在四倍量级到了四维全状态n 翻倍是十六倍。这不是线性关系选参数之前先算一遍内存再动手。参数物理含义常用取值调大的代价n每维网格点数101 / 201 / 401内存约 n^d四维时翻倍即 16 倍L计算域半边长初始距离的 1.5–2 倍同 n 下 dx 变粗界面精度下降cfl时间步系数0.3–0.5超过 0.5 一阶格式开始振荡alpha数值粘性系数取 |a−b|零水平集被抹厚捕获判定偏保守T逆向推进时长由业务问题定越长越容易吃到边界污染实操顺序我建议反过来先固定 L 和 T用 n51、101、201 各跑一遍看零水平集半径的误差随 dx 怎么降。一阶格式的误差应该大致按 dx 线性下降如果 101 到 201 误差几乎不动说明误差来源不是离散精度而是别的东西——多半是边界或者时间没走满回去查第 3.3 节那两个点。4.2 数值粘性不是越小越好很多人第一反应是“粘性就是误差能不能调到 0”。在这类双曲型方程上不行。粘性系数低于 CFL 允许的下限零水平集上会长出一圈高频振荡画出来像毛刺策略提取时梯度符号在相邻网格间来回翻闭环仿真直接发散。粘性偏大则是另一种难受界面被抹成一条宽带零水平集的位置往保守方向偏捕获区算小了你会得出“明明追得上却判定追不上”的结论。判断标准很朴素——拿简单运动的解析解对拍半径误差控制在两三个网格以内就算合格超了就先怀疑粘性系数是不是随手填的常数而不是跟着 |a−b| 走。还有个细节a 和 b 接近时α |a−b| 趋近于零格式本身退化。这不是 bug物理上确实如此——速度相等时捕获集不扩张。但数值上你会看到零水平集几乎不动、同时残留一点振荡这时候应该把问题重新表述比如加上时间上限改成“能否在 T 内追上”而不是硬把 α 调大凑出一个好看的结果。4.3 收敛判据与停止条件定长时间推进不需要收敛判据跑满 steps 就停。但如果你想求稳态无限时域的价值函数就得设停止条件常见做法是监视相邻两步的 L∞ 变化def solve_steady(n201, L2.0, a1.0, b0.6, cfl0.4, tol1e-4, max_steps20000): X, Y, dx build_grid(n, L) V terminal_value(X, Y) speed max(abs(a - b), 1e-12) dt cfl * dx / speed for k in range(max_steps): dV dt * hji_rhs(V, dx, a, b) V V dV V[0, :] V[1, :] dx; V[-1, :] V[-2, :] dx V[:, 0] V[:, 1] dx; V[:, -1] V[:, -2] dx # 用增量而不是整场变化量做判据, 避免大数值掩盖小区域的漂移 if np.max(np.abs(dV)) tol: break return X, Y, Vtol取 1e-4 配合 dx 在 0.01 量级是合理的比例tol 设得比 dx 还小两三个数量级没有意义只会白烧算力。另外注意判据用的是单步增量而不是整场的相对变化后者会被大数值区域拉平导致边界附近的缓慢漂移检测不到。5. 避坑与排查五个把人熬到后半夜的翻车现场5.1 边界零填充凭空造出一个方形捕获区现象零水平集本该是个圆结果跑出来贴着计算域边缘长成方框四个角还被拉出尖刺。原因边界值被固定成 0 或者某个常数等价于告诉求解器“边界处的价值是零”也就是“边界上就能完成任务”。捕获区会一路往外长直到被物理边界挡住。解决外边界的值必须按外推给V[边界] V[内邻] dx。同时做一个自检——跑完之后看零水平集有没有碰到计算域边界碰上了就把 L 放大或者把 T 缩短两个办法选一个。5.2 时间方向搞反捕获集反而在缩小现象a b 明明追捕方更快跑出来的零水平集却随时间收缩。原因H 的符号或者 min/max 的顺序写反了。H (b−a)|∇V| 里 b−a 是负的逆向推进方程是 V_τ −H两个负号叠起来才是扩张。中间任何一处翻符号结果就是反的。解决拿最简单的情形做单元测试——a1、b0.6、初始半径 0.2跑 τ 之后零水平集半径应该是 0.2 0.4τ误差在几个网格内。这个测试五秒钟就能写完能挡住后面 80% 的符号问题。5.3 量纲混着填追捕时间差了三倍多现象换了真实平台参数之后算出来的捕获时间跟实测对不上差一个固定倍数。原因一侧填的是 km/h另一侧填的是 m/sa/b 的比值被悄悄改掉了。这类问题在纯数值实验里不会暴露一接真实数据就现形。解决全部归一化。长度单位取初始相对距离时间单位取初始距离除以最大速度两个速度变成无量纲数你只需要关心 a/b 这个比值。归一化还有个附带好处网格参数不用每次跟着量纲重调。5.4 梯度噪声让策略在追与不追之间反复切换现象把反馈策略接到闭环仿真里追捕器在某个位置附近打转或者来回摆头。原因一个是前面说的 |∇V| 趋零没做保护另一个是中心差分在界面附近本身就有振荡梯度方向在相邻网格之间翻转。解决加 eps 下限只是止血根治要靠格式把一阶差分成对梯度做单调化处理或者直接上更高阶的加权重构再用双线性插值把网格上的策略插到连续状态去用。插值这一步别省直接用最近邻取策略闭环里会出现台阶式的抖振。5.5 四维网格直接把内存打爆现象换成全状态追捕方二维 逃逸方二维之后脚本跑到一半 MemoryError或者机器开始疯狂换页。原因四维网格点数按 n⁴ 长。n101 就是一亿多个点双精度 800MB 起加上中间量梯度、拉普拉斯、时间推进的临时数组实际占用要乘三四倍。解决优先做状态约化。简单运动下相对位置二维就够信息一点没丢确实需要各向异性转弯半径、加速度约束时先降 n 到 51 探路确认结果量级对了再往上加或者把数组换成单精度、把循环丢到 GPU 上。中间量用out参数复用同一块内存别每次都在函数里新建np.zeros_like。6. 验证与进阶拿解析解给数值解做校验数值解最怕“看着像那么回事”。简单运动这个算例的好处是有闭式解零水平集应该是一个半径 (a−b)τ 的圆价值函数就是 |x| − (a−b)τ。我一般做三组对照。第一组看半径。跑 τ 0.25、0.5、1.0 三个时刻量零水平集的最大半径跟解析值比。一阶格式的误差大致按 dx 线性下降量级参考网格 ndx半径相对误差估计收敛阶510.086% 上下—1010.043% 上下约 12010.021.5% 上下约 1如果三条误差不按这个趋势走别怀疑解析解回去查边界和粘性。第二组看闭环。用第 3.4 节提取的策略做一次前向积分初始距离 d捕获时间应该收敛到 d/(a−b)。这一组能抓出策略方向反了、插值没做、eps 太小这三类问题。第三组换个 a b 的算例。此时捕获集应该收缩初始不在集合内的点永远逃掉。很多符号错误在第一组里看不出来因为对称到这一组立刻暴露。进阶方向我建议从“加约束”入手给追捕方加转弯半径上限哈密顿量不再各向同性min over u 不能再解析地化简成 −a|∇V|得在离散层面枚举控制或者预计算每个方向的极值。这一步做完简单运动那套参考实现就不能用了但前两组的验证框架可以原样搬过去只是解析解换成“已知的屏障面形状”做定性比对。血泪经验是每次改动动力学项之后哪怕只改了一个系数也要先把简单运动的解析对照重跑一遍。有回偷懒跳过了符号错了半页闭环仿真里追捕器优雅地绕着目标转了三圈才被我看出不对。从那以后我每次动 HJI 那几行都强制走一遍半径对照。希望帮到你。本文还有配套的精品资源点击获取
