简介微分游戏difgames开源项目面向博弈论与控制理论方向的学习者和研究人员聚焦以哈密顿-雅可比-贝尔曼-伊萨克斯HJBI方程求解连续时间动态博弈中的纳什均衡问题适合具备一定数学功底与C语言基础、希望从代码层面理解动态优化与策略选择的读者。压缩包共9个文件、约206KB以6个C源文件为主体分别对应不同算法与仿真场景下的偏微分方程数值求解1个头文件负责绘制博弈过程与结果的图形界面1个Shell脚本用于自动化编译构建另有1份PDF文档阐述优化思路、算法原理与实验结论结构紧凑在Linux环境下即可编译运行。目前已有239人学习下载。读者可从中获得一套可运行的HJBI求解原型代码理解有限差分等数值近似手段如何落地并借助可视化输出直观验证策略收敛过程文档与源码相互印证也便于将动态博弈模型迁移到自动驾驶、经济建模等场景做二次实验。1. 微分游戏听起来像论文黑话但这套开源软件能回答“追不追得上、多久追上”微分游戏differential game中文里也常写微分博弈在开源软件里属于冷门品类论文一大堆能直接跑起来的工程没几个。我第一次接触它是因为一个追逃场景拦截方速度略快逃跑方机动性更好规则写出来不到十行但“到底能不能追上”“追上要几秒”“双方各自的最优动作长什么样”这三个问题靠调参和试凑永远回答不了。difgames 这套代码做的就是这件事把连续时间的对抗过程写成状态方程和控制集再用数值方法把值函数解出来最后从值函数的梯度反推双方的最优控制律。它适合做拦截与规避仿真、对抗策略基准、博弈论课程的动手部分也适合想给强化学习找一个“有解析解可对拍”的验证环境的人。2. 建模落笔顺序状态量、控制集、代价函数写错一个后面全白算建模这一步最容易出错因为它不报错。代码能跑、图也能出但值函数是错的而且错得很像对的。所以我把这套包的建模层拆成三个必须按顺序确认的东西状态量选什么坐标系、控制集是什么形状、终端代价怎么给。顺序不能颠倒因为后面每一步的数值方案都被前面定死了。2.1 状态量怎么选相对坐标能砍掉一半维度最直觉的写法是把追方和逃方位置都放进状态向量二维平面就是四维状态 $(p_x, p_y, e_x, e_y)$。我不建议一上来就这么干。如果动力学是平移不变的——也就是双方的速度上限不随绝对位置变化、场上没有障碍物和边界——那么值函数只依赖相对位置 $r p_{evader} - p_{pursuer}$四维问题直接塌缩成二维。这个判断有明确的适用条件不是偷懒控制集必须是常数速度圆盘不随位置变形、动力学里不能出现绝对坐标项比如随位置变化的风场、地形阻力、场上不能有障碍物或不可通行区。三条里破任意一条相对坐标就不成立必须回到全状态空间或者把破坏平移不变性的那部分额外塞进状态向量。我在第三章会讲这个代价有多大。判定标准很简单把相对动力学写出来如果 $\dot r$ 的表达式中不显含 $p$ 和 $e$相对坐标就是合法的。写出来显含了还硬用属于自欺欺人。2.2 动力学、控制集、代价函数的三段式写法这套包的约定是把一个对抗问题写成三段动力学 $\dot r f(r, u, v)$、控制约束 $u \in \mathcal{U}, v \in \mathcal{V}$、终端代价 $\phi(r)$。追方 $u$ 取最小、逃方 $v$ 取最大。本算例取最简单的全向运动模型$$\dot r v - u,\quad |u| \le a,\quad |v| \le b,\quad a b$$三段都落地之后哈密顿量可以手推出来$$H(\nabla V) \min_{|u|\le a}\max_{|v|\le b}\nabla V^\top (v-u) (b-a)|\nabla V|$$记速度优势 $s a - b$。$s 0$ 时捕获才有可能$s 0$ 时值函数根本不演化只有初始就落在捕获集里的状态才算捕获$s 0$ 时在无界域内逃方永远跑得掉。这个结论看起来平凡但我见过不少人第一版代码里把 a 和 b 的位置写反然后花了半天怀疑求解器。注意$s$ 的符号决定了整个方程的迎风方向。$s$ 恒正时求解器可以一直用同一套单侧差分如果场上存在 $s$ 变号的区域迎风方向要跟着符号翻转否则格式立刻失稳。2.3 最小可跑闭环定义 → 离散 → 迭代 → 出图下面的代码是建模层的完整落笔包含网格、控制集、终端代价三部分。参数我按“笔记本上几秒钟跑完”的量级设的改之前先把这一版跑通。import numpy as np # ---- 1. 网格状态是相对位置 r p_evader - p_pursuer覆盖 [-L, L]^2 ---- N 201 # 每维节点数取奇数保证原点落在节点上 L 5.0 # 相对位置的有效范围单位 m x np.linspace(-L, L, N) y np.linspace(-L, L, N) X, Y np.meshgrid(x, y, indexingij) dx, dy x[1] - x[0], y[1] - y[0] # 0.05 m # ---- 2. 控制集全向速度圆盘追方上限 a逃方上限 b ---- a, b 1.2, 0.8 # m/s s a - b # 速度优势必须 0 R_cap 0.3 # 捕获半径|r| R_cap 视为捕获成功 # ---- 3. 终端代价距离型初值捕获集内恒为 0 ---- # 不要写成 0/1 硬跳变边界处的无穷大梯度会把格式打爆 phi np.maximum(0.0, np.hypot(X, Y) - R_cap) capture_mask (np.hypot(X, Y) R_cap)逐段说明。第 1 段N201配L5得到dx0.05这个分辨率对R_cap0.3的捕获圆来说每半径有 6 个节点勉强够用再粗就会出现“捕获集在网格上是个十字形”的失真。第 2 段ab必须显式断言别赖在心里。第 3 段phi用距离而不是 0/1是因为初值的梯度决定了第一次迭代的推进方向硬跳变会让单侧差分算出一个巨大的值第一步就把解甩飞。参数含义调整时的约束N每维节点数与L共同决定dxdx应远小于R_capL相对位置范围至少要覆盖感兴趣的最远初始距离a/b双方速度上限a-b决定推进速度也决定 CFL 步长R_cap捕获半径小于3*dx时结果基本不可信跑通这一步之后应该能手算验证一个初始相对距离 $|r_0|$ 的状态捕获时间理论上是 $(|r_0| - R_{cap})/s$。这就是后面所有验证的基准线。3. HJI 方程的数值落地迎风差分、CFL 与网格无关性建模对了下一步才是数值。这套包的求解层是我最看重的部分因为它没有用“调包解 PDE”那种省事写法而是老老实实做了迎风格式和时间步长控制。微分博弈的值函数满足的是 Hamilton-Jacobi-Isaacs 方程一阶、非线性、解可能不光滑通用的 PDE 求解器在这里基本没用。3.1 为什么不去直接解那个偏微分方程HJI 方程写成演化形式是 $\partial_t V H(\nabla V) 0$配上终端条件 $V(r,0) \phi(r)$ 反向推进。麻烦在于 $H$ 只依赖 $\nabla V$方程是纯双曲的解在特征线相交处会形成激波——也就是值函数的梯度出现间断。在那里中心差分会振荡、二阶格式会过冲只有迎风类格式能收敛到正确的粘性解。另一条路是 fast marching / fast sweeping只解稳态的 eikonal 方程 $|\nabla T| 1/s$速度快得多。我一开始也想走这条捷径但它要求 $s$ 是常数一旦速度优势随状态变化第三章末和第六章都会碰到eikonal 形式就失效了而演化形式的格式几乎不用改。所以这套包选了演化形式作为主路径把快扫法留作特殊情形下的加速手段。3.2 迎风差分算对 |∇V| 才有后面的事核心就一个函数在网格上算迎风意义下的梯度模长。Godunov 格式的写法是每维取“后向差分为正”和“前向差分为负”里更大的那个再平方求和。def godunov_grad(V, dx, dy): 迎风模长 |∇V|用于 ∂V/∂t s|∇V| Dxm np.zeros_like(V); Dxm[1:, :] (V[1:, :] - V[:-1, :]) / dx # 后向 Dxp np.zeros_like(V); Dxp[:-1, :] (V[1:, :] - V[:-1, :]) / dx # 前向 Dym np.zeros_like(V); Dym[:, 1:] (V[:, 1:] - V[:, :-1]) / dy Dyp np.zeros_like(V); Dyp[:, :-1] (V[:, 1:] - V[:, :-1]) / dy # Godunov每个方向只取“会把信息送进来”的那一侧 gx np.maximum(np.maximum(Dxm, 0.0) ** 2, np.minimum(Dxp, 0.0) ** 2) gy np.maximum(np.maximum(Dym, 0.0) ** 2, np.minimum(Dyp, 0.0) ** 2) return np.sqrt(gx gy)Dxm是后向差分、Dxp是前向差分边界上它们的初值保持 0这样域外不会引入鬼影数据。np.maximum(Dxm, 0)表示只有当后向差分为正值函数在向左增长时左侧的信息才是有效来源np.minimum(Dxp, 0)对应右侧。两个方向平方取大再相加正是 Godunov 格式对 $|\nabla V|$ 的标准离散。一个常见误解是“用中心差分更精确”。在这里恰恰相反中心差分在光滑区确实二阶但在激波处会产生 $O(1)$ 的过冲而且过冲会被时间推进不断放大最后表现为值函数表面出现细密的高频纹路。这套包的策略是全场统一用迎风牺牲光滑区的精度换稳定性——值函数这东西稳定比精确重要得多。3.3 步长与 CFLdt 取大了会看到棋盘格显式推进的稳定性由 CFL 条件约束。对多维迎风格式粗略的约束是 $dt \cdot s \cdot (1/dx 1/dy) \le 1$工程上再乘一个 0.40.5 的安全系数。s_max s # 本算例 s 恒定变速情形取全场最大值 dt 0.4 * min(dx, dy) / s_max # 0.4*0.05/0.4 0.05 s n_iter 2000 # 波前穿过半场约需 L/s/dt 250 步留足余量 V phi.copy() for k in range(n_iter): g godunov_grad(V, dx, dy) Vn V dt * s * g Vn[capture_mask] 0.0 # 捕获集内部钉死为 0防止数值扩散渗进来 V Vn if k % 200 0: t k * dt exact np.maximum(0.0, np.hypot(X, Y) - R_cap - s * t) print(fiter{k:4d} max|V-exact|{np.abs(V-exact).max():.2e})dt那一行的逻辑min(dx,dy)取最不利的网格间距除以s_max保证单步推进不超过一个网格再乘 0.4 留余量。n_iter的估算用L/s/dt——波前从中心跑到边界需要的步数取十倍左右比较保险。Vn[capture_mask] 0.0这一行看着多余其实不能少。数值扩散会让捕获集内部的 $V$ 慢慢变成一个小正数累积几百步之后零等值线会往里缩捕获区就凭空变小了。最后那段打印是本算例的“后悔药”因为 $s$ 恒定、初值是距离函数解析解是 $V(r,t) \max(0, |r| - R_{cap} - st)$可以直接对拍。残差应该在一两个网格步长对应的量级大约 $10^{-2}$如果第 200 步就跑到 $10^{-1}$ 以上说明步长或格式有问题别再往下调参数了先回来查这两处。提示波前一旦触到网格边界解就不可信了。加一句print(V[0,:].max(), V[-1,:].max())盯着边界上的值是不是还在单调增长比看整张图更早发现问题。4. 从值函数到控制律策略提取与 200 步闭环回放验证值函数解出来只是中间产物真正要用的是策略。这一步在微分游戏里比在最优控制里更需要小心因为双方策略都从同一个 $\nabla V$ 里取符号和方向错一个整个验证就变成自娱自乐。4.1 梯度方向决定最优控制对 §2.2 的哈密顿量做 argmin / argmax$$u^* a\frac{\nabla V}{|\nabla V|},\qquad v^* b\frac{\nabla V}{|\nabla V|}$$两者同向这看起来反直觉但在相对坐标下是对的$\dot r v - u$追方沿 $\nabla V$ 走贡献是 $-a\hat g$把 $r$ 往捕获集里推逃方沿同方向走贡献 $b\hat g$往外拉。净效果是 $r$ 以 $s$ 的速率沿 $-\hat g$ 收缩。换成绝对坐标写追方的运动方向就变成 $-\nabla V$ 了——这个符号翻转是踩坑重灾区第五章会单独说。实现上梯度必须和哈密顿量用同一套迎风逻辑不能换成中心差分def grad_upwind(V, dx, dy): 与哈密顿量一致的梯度向量每维取模长更大的单侧差分 Dxm np.zeros_like(V); Dxm[1:, :] (V[1:, :] - V[:-1, :]) / dx Dxp np.zeros_like(V); Dxp[:-1, :] (V[1:, :] - V[:-1, :]) / dx Dym np.zeros_like(V); Dym[:, 1:] (V[:, 1:] - V[:, :-1]) / dy Dyp np.zeros_like(V); Dyp[:, :-1] (V[:, 1:] - V[:, :-1]) / dy Vx np.where(np.abs(Dxm) np.abs(Dxp), Dxm, Dxp) Vy np.where(np.abs(Dym) np.abs(Dyp), Dym, Dyp) return Vx, Vynp.where那一行是“取模长更大的一侧”等价于 Godunov 选出的那个方向只是这里需要保留符号所以不能像算模长时那样直接取平方。梯度场Vx, Vy只算一次并缓存闭环仿真里反复插值取用。4.2 闭环回放200 步仿真看三个指标策略有了把双方按各自的最优控制律同时施加在相对坐标里积分就能回放整场对抗。这套包的验证层就是干这个的。def interp2(Vx, Vy, x, y, p): 双线性插值取梯度p(px,py) 为相对位置 i int(np.clip(np.searchsorted(x, p[0]) - 1, 0, len(x) - 2)) j int(np.clip(np.searchsorted(y, p[1]) - 1, 0, len(y) - 2)) fx (p[0] - x[i]) / (x[1] - x[0]); fy (p[1] - y[j]) / (y[1] - y[0]) w np.array([(1-fx)*(1-fy), fx*(1-fy), (1-fx)*fy, fx*fy]) gx w np.array([Vx[i, j], Vx[i1, j], Vx[i, j1], Vx[i1, j1]]) gy w np.array([Vy[i, j], Vy[i1, j], Vy[i, j1], Vy[i1, j1]]) return np.array([gx, gy]) def rollout(Vx, Vy, x, y, r0, a, b, dt0.01, tmax30.0): r np.array(r0, dtypefloat); path [r.copy()] for k in range(int(tmax / dt)): g interp2(Vx, Vy, x, y, r) n np.linalg.norm(g) if n 1e-9: # 落进值函数平底区梯度消失 break u a * g / n v b * g / n r r dt * (v - u) # 相对坐标净收缩速度 s path.append(r.copy()) if np.linalg.norm(r) R_cap: return k * dt, True, np.array(path) return tmax, False, np.array(path)dt0.01是积分步长和求解器的dt0.05不是一回事可以更细tmax是放弃阈值超过就判定逃逸成功。返回值三元组是捕获时间、是否捕获、轨迹。三个要看的指标捕获时间和 $(|r_0|-R_{cap})/s$ 的偏差应该在网格分辨率量级、轨迹是否单调逼近原点出现来回横跳说明梯度有噪声、控制量 $|u|$ 是否稳定贴着 $a$频繁在 0 和 a 之间跳变就是抖振。4.3 网格无关性三组分辨率下结果该收敛到什么程度单次跑通说明不了什么我习惯用三组分辨率交叉验证。以 $r_0(3,2)$ 为例$|r_0|3.606$理论捕获时间 $(3.606-0.3)/0.4 8.26$ 秒。网格dx仿真捕获时间相对误差单次迭代耗时101×1010.108.38 s本文还有配套的精品资源点击获取 
长河君 · 2026/10/29 14:23:16 · 阅读 1.2k 微分博弈HJI数值求解：追逃策略提取与闭环验证
简介微分游戏difgames开源项目面向博弈论与控制理论方向的学习者和研究人员聚焦以哈密顿-雅可比-贝尔曼-伊萨克斯HJBI方程求解连续时间动态博弈中的纳什均衡问题适合具备一定数学功底与C语言基础、希望从代码层面理解动态优化与策略选择的读者。压缩包共9个文件、约206KB以6个C源文件为主体分别对应不同算法与仿真场景下的偏微分方程数值求解1个头文件负责绘制博弈过程与结果的图形界面1个Shell脚本用于自动化编译构建另有1份PDF文档阐述优化思路、算法原理与实验结论结构紧凑在Linux环境下即可编译运行。目前已有239人学习下载。读者可从中获得一套可运行的HJBI求解原型代码理解有限差分等数值近似手段如何落地并借助可视化输出直观验证策略收敛过程文档与源码相互印证也便于将动态博弈模型迁移到自动驾驶、经济建模等场景做二次实验。1. 微分游戏听起来像论文黑话但这套开源软件能回答“追不追得上、多久追上”微分游戏differential game中文里也常写微分博弈在开源软件里属于冷门品类论文一大堆能直接跑起来的工程没几个。我第一次接触它是因为一个追逃场景拦截方速度略快逃跑方机动性更好规则写出来不到十行但“到底能不能追上”“追上要几秒”“双方各自的最优动作长什么样”这三个问题靠调参和试凑永远回答不了。difgames 这套代码做的就是这件事把连续时间的对抗过程写成状态方程和控制集再用数值方法把值函数解出来最后从值函数的梯度反推双方的最优控制律。它适合做拦截与规避仿真、对抗策略基准、博弈论课程的动手部分也适合想给强化学习找一个“有解析解可对拍”的验证环境的人。2. 建模落笔顺序状态量、控制集、代价函数写错一个后面全白算建模这一步最容易出错因为它不报错。代码能跑、图也能出但值函数是错的而且错得很像对的。所以我把这套包的建模层拆成三个必须按顺序确认的东西状态量选什么坐标系、控制集是什么形状、终端代价怎么给。顺序不能颠倒因为后面每一步的数值方案都被前面定死了。2.1 状态量怎么选相对坐标能砍掉一半维度最直觉的写法是把追方和逃方位置都放进状态向量二维平面就是四维状态 $(p_x, p_y, e_x, e_y)$。我不建议一上来就这么干。如果动力学是平移不变的——也就是双方的速度上限不随绝对位置变化、场上没有障碍物和边界——那么值函数只依赖相对位置 $r p_{evader} - p_{pursuer}$四维问题直接塌缩成二维。这个判断有明确的适用条件不是偷懒控制集必须是常数速度圆盘不随位置变形、动力学里不能出现绝对坐标项比如随位置变化的风场、地形阻力、场上不能有障碍物或不可通行区。三条里破任意一条相对坐标就不成立必须回到全状态空间或者把破坏平移不变性的那部分额外塞进状态向量。我在第三章会讲这个代价有多大。判定标准很简单把相对动力学写出来如果 $\dot r$ 的表达式中不显含 $p$ 和 $e$相对坐标就是合法的。写出来显含了还硬用属于自欺欺人。2.2 动力学、控制集、代价函数的三段式写法这套包的约定是把一个对抗问题写成三段动力学 $\dot r f(r, u, v)$、控制约束 $u \in \mathcal{U}, v \in \mathcal{V}$、终端代价 $\phi(r)$。追方 $u$ 取最小、逃方 $v$ 取最大。本算例取最简单的全向运动模型$$\dot r v - u,\quad |u| \le a,\quad |v| \le b,\quad a b$$三段都落地之后哈密顿量可以手推出来$$H(\nabla V) \min_{|u|\le a}\max_{|v|\le b}\nabla V^\top (v-u) (b-a)|\nabla V|$$记速度优势 $s a - b$。$s 0$ 时捕获才有可能$s 0$ 时值函数根本不演化只有初始就落在捕获集里的状态才算捕获$s 0$ 时在无界域内逃方永远跑得掉。这个结论看起来平凡但我见过不少人第一版代码里把 a 和 b 的位置写反然后花了半天怀疑求解器。注意$s$ 的符号决定了整个方程的迎风方向。$s$ 恒正时求解器可以一直用同一套单侧差分如果场上存在 $s$ 变号的区域迎风方向要跟着符号翻转否则格式立刻失稳。2.3 最小可跑闭环定义 → 离散 → 迭代 → 出图下面的代码是建模层的完整落笔包含网格、控制集、终端代价三部分。参数我按“笔记本上几秒钟跑完”的量级设的改之前先把这一版跑通。import numpy as np # ---- 1. 网格状态是相对位置 r p_evader - p_pursuer覆盖 [-L, L]^2 ---- N 201 # 每维节点数取奇数保证原点落在节点上 L 5.0 # 相对位置的有效范围单位 m x np.linspace(-L, L, N) y np.linspace(-L, L, N) X, Y np.meshgrid(x, y, indexingij) dx, dy x[1] - x[0], y[1] - y[0] # 0.05 m # ---- 2. 控制集全向速度圆盘追方上限 a逃方上限 b ---- a, b 1.2, 0.8 # m/s s a - b # 速度优势必须 0 R_cap 0.3 # 捕获半径|r| R_cap 视为捕获成功 # ---- 3. 终端代价距离型初值捕获集内恒为 0 ---- # 不要写成 0/1 硬跳变边界处的无穷大梯度会把格式打爆 phi np.maximum(0.0, np.hypot(X, Y) - R_cap) capture_mask (np.hypot(X, Y) R_cap)逐段说明。第 1 段N201配L5得到dx0.05这个分辨率对R_cap0.3的捕获圆来说每半径有 6 个节点勉强够用再粗就会出现“捕获集在网格上是个十字形”的失真。第 2 段ab必须显式断言别赖在心里。第 3 段phi用距离而不是 0/1是因为初值的梯度决定了第一次迭代的推进方向硬跳变会让单侧差分算出一个巨大的值第一步就把解甩飞。参数含义调整时的约束N每维节点数与L共同决定dxdx应远小于R_capL相对位置范围至少要覆盖感兴趣的最远初始距离a/b双方速度上限a-b决定推进速度也决定 CFL 步长R_cap捕获半径小于3*dx时结果基本不可信跑通这一步之后应该能手算验证一个初始相对距离 $|r_0|$ 的状态捕获时间理论上是 $(|r_0| - R_{cap})/s$。这就是后面所有验证的基准线。3. HJI 方程的数值落地迎风差分、CFL 与网格无关性建模对了下一步才是数值。这套包的求解层是我最看重的部分因为它没有用“调包解 PDE”那种省事写法而是老老实实做了迎风格式和时间步长控制。微分博弈的值函数满足的是 Hamilton-Jacobi-Isaacs 方程一阶、非线性、解可能不光滑通用的 PDE 求解器在这里基本没用。3.1 为什么不去直接解那个偏微分方程HJI 方程写成演化形式是 $\partial_t V H(\nabla V) 0$配上终端条件 $V(r,0) \phi(r)$ 反向推进。麻烦在于 $H$ 只依赖 $\nabla V$方程是纯双曲的解在特征线相交处会形成激波——也就是值函数的梯度出现间断。在那里中心差分会振荡、二阶格式会过冲只有迎风类格式能收敛到正确的粘性解。另一条路是 fast marching / fast sweeping只解稳态的 eikonal 方程 $|\nabla T| 1/s$速度快得多。我一开始也想走这条捷径但它要求 $s$ 是常数一旦速度优势随状态变化第三章末和第六章都会碰到eikonal 形式就失效了而演化形式的格式几乎不用改。所以这套包选了演化形式作为主路径把快扫法留作特殊情形下的加速手段。3.2 迎风差分算对 |∇V| 才有后面的事核心就一个函数在网格上算迎风意义下的梯度模长。Godunov 格式的写法是每维取“后向差分为正”和“前向差分为负”里更大的那个再平方求和。def godunov_grad(V, dx, dy): 迎风模长 |∇V|用于 ∂V/∂t s|∇V| Dxm np.zeros_like(V); Dxm[1:, :] (V[1:, :] - V[:-1, :]) / dx # 后向 Dxp np.zeros_like(V); Dxp[:-1, :] (V[1:, :] - V[:-1, :]) / dx # 前向 Dym np.zeros_like(V); Dym[:, 1:] (V[:, 1:] - V[:, :-1]) / dy Dyp np.zeros_like(V); Dyp[:, :-1] (V[:, 1:] - V[:, :-1]) / dy # Godunov每个方向只取“会把信息送进来”的那一侧 gx np.maximum(np.maximum(Dxm, 0.0) ** 2, np.minimum(Dxp, 0.0) ** 2) gy np.maximum(np.maximum(Dym, 0.0) ** 2, np.minimum(Dyp, 0.0) ** 2) return np.sqrt(gx gy)Dxm是后向差分、Dxp是前向差分边界上它们的初值保持 0这样域外不会引入鬼影数据。np.maximum(Dxm, 0)表示只有当后向差分为正值函数在向左增长时左侧的信息才是有效来源np.minimum(Dxp, 0)对应右侧。两个方向平方取大再相加正是 Godunov 格式对 $|\nabla V|$ 的标准离散。一个常见误解是“用中心差分更精确”。在这里恰恰相反中心差分在光滑区确实二阶但在激波处会产生 $O(1)$ 的过冲而且过冲会被时间推进不断放大最后表现为值函数表面出现细密的高频纹路。这套包的策略是全场统一用迎风牺牲光滑区的精度换稳定性——值函数这东西稳定比精确重要得多。3.3 步长与 CFLdt 取大了会看到棋盘格显式推进的稳定性由 CFL 条件约束。对多维迎风格式粗略的约束是 $dt \cdot s \cdot (1/dx 1/dy) \le 1$工程上再乘一个 0.40.5 的安全系数。s_max s # 本算例 s 恒定变速情形取全场最大值 dt 0.4 * min(dx, dy) / s_max # 0.4*0.05/0.4 0.05 s n_iter 2000 # 波前穿过半场约需 L/s/dt 250 步留足余量 V phi.copy() for k in range(n_iter): g godunov_grad(V, dx, dy) Vn V dt * s * g Vn[capture_mask] 0.0 # 捕获集内部钉死为 0防止数值扩散渗进来 V Vn if k % 200 0: t k * dt exact np.maximum(0.0, np.hypot(X, Y) - R_cap - s * t) print(fiter{k:4d} max|V-exact|{np.abs(V-exact).max():.2e})dt那一行的逻辑min(dx,dy)取最不利的网格间距除以s_max保证单步推进不超过一个网格再乘 0.4 留余量。n_iter的估算用L/s/dt——波前从中心跑到边界需要的步数取十倍左右比较保险。Vn[capture_mask] 0.0这一行看着多余其实不能少。数值扩散会让捕获集内部的 $V$ 慢慢变成一个小正数累积几百步之后零等值线会往里缩捕获区就凭空变小了。最后那段打印是本算例的“后悔药”因为 $s$ 恒定、初值是距离函数解析解是 $V(r,t) \max(0, |r| - R_{cap} - st)$可以直接对拍。残差应该在一两个网格步长对应的量级大约 $10^{-2}$如果第 200 步就跑到 $10^{-1}$ 以上说明步长或格式有问题别再往下调参数了先回来查这两处。提示波前一旦触到网格边界解就不可信了。加一句print(V[0,:].max(), V[-1,:].max())盯着边界上的值是不是还在单调增长比看整张图更早发现问题。4. 从值函数到控制律策略提取与 200 步闭环回放验证值函数解出来只是中间产物真正要用的是策略。这一步在微分游戏里比在最优控制里更需要小心因为双方策略都从同一个 $\nabla V$ 里取符号和方向错一个整个验证就变成自娱自乐。4.1 梯度方向决定最优控制对 §2.2 的哈密顿量做 argmin / argmax$$u^* a\frac{\nabla V}{|\nabla V|},\qquad v^* b\frac{\nabla V}{|\nabla V|}$$两者同向这看起来反直觉但在相对坐标下是对的$\dot r v - u$追方沿 $\nabla V$ 走贡献是 $-a\hat g$把 $r$ 往捕获集里推逃方沿同方向走贡献 $b\hat g$往外拉。净效果是 $r$ 以 $s$ 的速率沿 $-\hat g$ 收缩。换成绝对坐标写追方的运动方向就变成 $-\nabla V$ 了——这个符号翻转是踩坑重灾区第五章会单独说。实现上梯度必须和哈密顿量用同一套迎风逻辑不能换成中心差分def grad_upwind(V, dx, dy): 与哈密顿量一致的梯度向量每维取模长更大的单侧差分 Dxm np.zeros_like(V); Dxm[1:, :] (V[1:, :] - V[:-1, :]) / dx Dxp np.zeros_like(V); Dxp[:-1, :] (V[1:, :] - V[:-1, :]) / dx Dym np.zeros_like(V); Dym[:, 1:] (V[:, 1:] - V[:, :-1]) / dy Dyp np.zeros_like(V); Dyp[:, :-1] (V[:, 1:] - V[:, :-1]) / dy Vx np.where(np.abs(Dxm) np.abs(Dxp), Dxm, Dxp) Vy np.where(np.abs(Dym) np.abs(Dyp), Dym, Dyp) return Vx, Vynp.where那一行是“取模长更大的一侧”等价于 Godunov 选出的那个方向只是这里需要保留符号所以不能像算模长时那样直接取平方。梯度场Vx, Vy只算一次并缓存闭环仿真里反复插值取用。4.2 闭环回放200 步仿真看三个指标策略有了把双方按各自的最优控制律同时施加在相对坐标里积分就能回放整场对抗。这套包的验证层就是干这个的。def interp2(Vx, Vy, x, y, p): 双线性插值取梯度p(px,py) 为相对位置 i int(np.clip(np.searchsorted(x, p[0]) - 1, 0, len(x) - 2)) j int(np.clip(np.searchsorted(y, p[1]) - 1, 0, len(y) - 2)) fx (p[0] - x[i]) / (x[1] - x[0]); fy (p[1] - y[j]) / (y[1] - y[0]) w np.array([(1-fx)*(1-fy), fx*(1-fy), (1-fx)*fy, fx*fy]) gx w np.array([Vx[i, j], Vx[i1, j], Vx[i, j1], Vx[i1, j1]]) gy w np.array([Vy[i, j], Vy[i1, j], Vy[i, j1], Vy[i1, j1]]) return np.array([gx, gy]) def rollout(Vx, Vy, x, y, r0, a, b, dt0.01, tmax30.0): r np.array(r0, dtypefloat); path [r.copy()] for k in range(int(tmax / dt)): g interp2(Vx, Vy, x, y, r) n np.linalg.norm(g) if n 1e-9: # 落进值函数平底区梯度消失 break u a * g / n v b * g / n r r dt * (v - u) # 相对坐标净收缩速度 s path.append(r.copy()) if np.linalg.norm(r) R_cap: return k * dt, True, np.array(path) return tmax, False, np.array(path)dt0.01是积分步长和求解器的dt0.05不是一回事可以更细tmax是放弃阈值超过就判定逃逸成功。返回值三元组是捕获时间、是否捕获、轨迹。三个要看的指标捕获时间和 $(|r_0|-R_{cap})/s$ 的偏差应该在网格分辨率量级、轨迹是否单调逼近原点出现来回横跳说明梯度有噪声、控制量 $|u|$ 是否稳定贴着 $a$频繁在 0 和 a 之间跳变就是抖振。4.3 网格无关性三组分辨率下结果该收敛到什么程度单次跑通说明不了什么我习惯用三组分辨率交叉验证。以 $r_0(3,2)$ 为例$|r_0|3.606$理论捕获时间 $(3.606-0.3)/0.4 8.26$ 秒。网格dx仿真捕获时间相对误差单次迭代耗时101×1010.108.38 s本文还有配套的精品资源点击获取