简介：微分游戏（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、谁在 max：HJI 的符号约定

两方零和，追捕方 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）上会出现折痕，最优策略在那里切换，用特征线法硬推会在折痕处直接崩掉，而网格法 + 数值粘性恰好收敛到粘性解，这是这套方法最值钱的地方。

代价是维数灾难。二维状态网格 n=201 是四万个点，跑起来很轻；一旦你把追捕方和逃逸方各自的二维位置都塞进网格，变成四维，n=81 就是四千三百万个点，双精度直接 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, indexing="ij") 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(n=201, L=2.0, a=1.0, b=0.6, T=1.5, cfl=0.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, eps=1e-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, vy

eps这个保护项是必需的。零水平集上 |∇V| 理论上不为零（距离函数的梯度恒为 1），但数值解在原点附近会因为网格对称性出现极小值，直接除会得到方向乱跳的单位向量，闭环仿真里表现为追捕器在目标附近抽搐。加了下限之后，靠近原点的方向精度下降，但至少不会给出反向指令。

跑完这三段，你手上应该有一张零水平集图：一个以 (a−b)τ 为半径向外扩张的圆。圆不圆、扩张速度对不对，就是下一章参数标定的起点。

4. 参数标定：dt、CFL 系数、数值粘性与收敛判据怎么定

4.1 网格分辨率与内存的取舍

先说清代价。二维网格下，n 翻倍内存是四倍、计算量也在四倍量级；到了四维全状态，n 翻倍是十六倍。这不是线性关系，选参数之前先算一遍内存再动手。

实操顺序我建议反过来：先固定 L 和 T，用 n=51、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(n=201, L=2.0, a=1.0, b=0.6, cfl=0.4, tol=1e-4, max_steps=20000): 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, V

tol取 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，两个负号叠起来才是扩张。中间任何一处翻符号，结果就是反的。

解决：拿最简单的情形做单元测试——a=1、b=0.6、初始半径 0.2，跑 τ 之后零水平集半径应该是 0.2 + 0.4τ，误差在几个网格内。这个测试五秒钟就能写完，能挡住后面 80% 的符号问题。

5.3 量纲混着填，追捕时间差了三倍多

现象：换了真实平台参数之后，算出来的捕获时间跟实测对不上，差一个固定倍数。

原因：一侧填的是 km/h，另一侧填的是 m/s，a/b 的比值被悄悄改掉了。这类问题在纯数值实验里不会暴露，一接真实数据就现形。

解决：全部归一化。长度单位取初始相对距离，时间单位取初始距离除以最大速度，两个速度变成无量纲数，你只需要关心 a/b 这个比值。归一化还有个附带好处：网格参数不用每次跟着量纲重调。

5.4 梯度噪声让策略在追与不追之间反复切换

现象：把反馈策略接到闭环仿真里，追捕器在某个位置附近打转，或者来回摆头。

原因：一个是前面说的 |∇V| 趋零没做保护，另一个是中心差分在界面附近本身就有振荡，梯度方向在相邻网格之间翻转。

解决：加 eps 下限只是止血；根治要靠格式，把一阶差分成对梯度做单调化处理（或者直接上更高阶的加权重构），再用双线性插值把网格上的策略插到连续状态去用。插值这一步别省，直接用最近邻取策略，闭环里会出现台阶式的抖振。

5.5 四维网格直接把内存打爆

现象：换成全状态（追捕方二维 + 逃逸方二维）之后，脚本跑到一半 MemoryError，或者机器开始疯狂换页。

原因：四维网格点数按 n⁴ 长。n=101 就是一亿多个点，双精度 800MB 起，加上中间量（梯度、拉普拉斯、时间推进的临时数组）实际占用要乘三四倍。

解决：优先做状态约化。简单运动下相对位置二维就够，信息一点没丢；确实需要各向异性（转弯半径、加速度约束）时，先降 n 到 51 探路，确认结果量级对了再往上加，或者把数组换成单精度、把循环丢到 GPU 上。中间量用out=参数复用同一块内存，别每次都在函数里新建np.zeros_like。

6. 验证与进阶：拿解析解给数值解做校验

数值解最怕“看着像那么回事”。简单运动这个算例的好处是有闭式解：零水平集应该是一个半径 (a−b)τ 的圆，价值函数就是 |x| − (a−b)τ。我一般做三组对照。

第一组看半径。跑 τ = 0.25、0.5、1.0 三个时刻，量零水平集的最大半径，跟解析值比。一阶格式的误差大致按 dx 线性下降，量级参考：

如果三条误差不按这个趋势走，别怀疑解析解，回去查边界和粘性。

第二组看闭环。用第 3.4 节提取的策略做一次前向积分，初始距离 d，捕获时间应该收敛到 d/(a−b)。这一组能抓出策略方向反了、插值没做、eps 太小这三类问题。

第三组换个 a < b 的算例。此时捕获集应该收缩，初始不在集合内的点永远逃掉。很多符号错误在第一组里看不出来（因为对称），到这一组立刻暴露。

进阶方向我建议从“加约束”入手：给追捕方加转弯半径上限，哈密顿量不再各向同性，min over u 不能再解析地化简成 −a|∇V|，得在离散层面枚举控制或者预计算每个方向的极值。这一步做完，简单运动那套参考实现就不能用了，但前两组的验证框架可以原样搬过去，只是解析解换成“已知的屏障面形状”做定性比对。

血泪经验是：每次改动动力学项之后，哪怕只改了一个系数，也要先把简单运动的解析对照重跑一遍。有回偷懒跳过了，符号错了半页，闭环仿真里追捕器优雅地绕着目标转了三圈才被我看出不对。从那以后我每次动 HJI 那几行，都强制走一遍半径对照。希望帮到你。

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