简介difgames 是一份面向博弈论、控制理论与科学计算学习者的开源代码包聚焦连续时间动态系统下的微分博弈问题借助哈密顿-雅可比-贝尔曼-伊萨克斯方程求解多决策者动态环境中的纳什均衡。资源共9个文件以6个C源文件为核心辅以1个头文件、1个Shell构建脚本和1份PDF说明文档压缩包约206KB源码涉及偏微分方程的数值近似与图形显示构建脚本可简化本地编译流程文档则补充优化方法与算法原理。目前已有240人学习。读者可据此研读HJBI方程的C语言实现路径理解有限差分等数值方法如何落地并通过图形界面观察博弈过程与结果进而把握动态博弈的建模、求解与验证思路相关框架也可为自动驾驶、经济模型和策略分析等场景提供参考。1. difgames把微分游戏从论文公式推到能跑的开源软件追逃、拦截、两车对抗这类问题在控制圈叫追踪-逃逸在博弈圈叫微分游戏——状态连续演化、双方同时下注、代价由整条轨迹决定。difgames 这个开源软件要解决的正是中间那段最难受的路公式推导完到可视化看结果之间得自己搭动力学、写数值格式、处理边界、抽策略、画轨迹。它把这些环节串成一条能跑通的链让一个二元对抗问题在几十分钟内出第一张值函数曲面。适合做对抗决策、多智能体博弈、控制课程设计的人也适合想验证某个解析结论到底对不对的算法岗。读完这篇重点是把微分游戏与开源软件两个字眼从概念落到代码。需要提前说清楚边界它不是游戏引擎做 NPC 的那种“游戏”也不是强化学习那套采样训练。它的核心是一条偏微分方程——值函数满足的 HJI 方程解出来才叫策略。所以能不能用取决于你的状态维度够不够低、动力学够不够光滑。2. 微分博弈的骨架四元组、相对坐标与三条求解路线拿到一个对抗问题先别急着写代码。我见过太多人上来就搭环境跑到一半发现动力学模型漏了速度约束值函数算出来是一张斜率离谱的斜面只能说重来。这一章把纸面上的东西先固定下来后面所有代码和参数的讨论都挂在这套骨架上。2.1 一个博弈在代码里必须落成哪四个对象微分博弈的标准写法是状态 $x\in\mathbb{R}^n$ 按 $\dot{x}f(t,x,u,v)$ 演化追方控制 $u\in\mathcal{U}$ 想最小化代价逃方控制 $v\in\mathcal{V}$ 想最大化代价代价泛函写成终端项加积分项$$J(u,v)\varphi(x(T))\int_0^T L(t,x,u,v),dt$$值函数 $V(t,x)\min_u\max_v J$ 满足的 HJI 方程是 $\partial_t V\min_u\max_v{\nabla V\cdot fL}0$终端条件 $V(T,x)\varphi(x)$。这一串符号落到工程上就是四个必须显式定义的对象缺一个都跑不起来。博弈要素数学对象代码里的落点写错的典型后果状态$x\in\mathbb{R}^n$环境的 state / 相对坐标维度和网格对不上广播错误不报错但结果全错动力学$\dot{x}f(t,x,u,v)$dynamics(t, x, u, v)漏加速度上限轨迹能瞬移代价$\varphi\int L$terminal_cost/running_cost量纲不匹配终端项被积分项压死控制集$|u|\le u_{max}$、$|v|\le v_{max}$环境里的速度上限$u_{max}\le v_{max}$ 时值函数无界求解器空转控制集那一行是最容易翻车的地方。追方速度上限只要不大于逃方在无界平面上就是永远追不上代码不会报错它只会让值函数在几万步迭代里线性地往上爬你以为是在收敛其实是在发散。我在第一次跑的时候就吃了这个亏日志里残差一直降不下去后来把 $u_{max}$ 从 0.8 调到 1.0三分钟就收敛了。2.2 相对坐标降维四维状态怎么压成两维网格追逃问题的自然状态是四个量追方位置两维、逃方位置两维。四维网格求解的代价是灾难性的161 个点每维就是 $161^4\approx6.7$ 亿个格点普通机器连一个时间步都扛不住。但这类问题几乎都有平移不变性把状态换成相对位置 $rp_{追}-p_{逃}$动力学变成 $\dot{r}u-v$状态只剩两维运行代价直接掉到 $161^2\approx2.6$ 万个格点单机几秒钟一步。这个变换不是免费的。它要求代价函数只依赖相对位置和相对速度一旦你的问题里有绝对位置相关的项比如边界约束、地形代价、某个固定目标点降维就不成立得回到全状态空间或者用别的分解方式。判断方法很朴素把动力学和代价里的绝对坐标全部替换成相对坐标如果表达式中不再出现单独的 $p_{追}$ 或 $p_{逃}$就可以降。降维之后还有一个隐含假设相对坐标平面是无限的。实际实现必须截断成一个有界方框然后在边界上做处理。这个截断大小直接决定你的求解域够不够覆盖实际轨迹取小了轨迹会贴边走取大了网格分辨率被稀释梯度估计噪声变大。我一般的做法是先用解析上界估一个粗的捕获距离再乘 2.5 倍当半边长。2.3 三条求解路线怎么选矩阵方程、网格 PDE、迭代最佳响应微分博弈的数值解法不是只有一条路选错路线比调错参数更浪费时间。下面这张表是我自己在不同项目里踩过之后总结的选型依据覆盖了绝大多数中小规模场景。路线适用条件精度计算成本主要风险线性二次型矩阵微分方程动力学线性、代价二次、无控制约束高近似解析极低秒级一旦有控制饱和就失效结论不可外推HJI 网格 PDE 推进状态维数 ≤ 3动力学光滑中受网格和格式限制中到高分钟级边界条件和迎风方向写错会出锯齿或伪解迭代最佳响应 / 打靶状态维数 4 以上只要开环轨迹低局部收敛低秒级收敛到局部均衡不保证全局鞍点线性二次型那条路值得单独提一句。它把 HJI 方程退化成一组矩阵微分方程从终端时间往回想推就能拿到时变反馈增益不需要任何网格也不需要处理差分格式。很多课程设计题、入门验证题都落在这个范围内能解析求解就别上网格纯属给自己找麻烦。需要网格 PDE 的场景一般有两个特征控制有饱和约束或者代价函数非二次比如带捕获半径的开关型终端代价。这时候 HJI 没有闭式解只能数值推进。推进的稳定性、精度、边界处理就是第 4 章要展开的内容。迭代最佳响应适合快速看个大概特别是状态维数到了 6 维以上网格法完全不可行的时候。它的思路是固定一方策略、优化另一方交替迭代到策略不再变化。代价是它给的是开环或局部反馈不是真正的值函数也没法画那张好看的曲面图。用它做前期探索可以出结论之前最好用低维网格解交叉验证一次。3. 把 difgames 跑起来从依赖自检到第一条追逃轨迹理论骨架立住了就可以动手。这一章按我实际的执行顺序走先确认环境、再定义博弈、然后解值函数、最后抽策略跑轨迹。每一步都有能直接抄的代码参数含义写在代码后面。3.1 依赖自检与最小可运行环境这类数值包的核心依赖其实很少数组运算、常微分方程求解、绘图。不需要深度学习框架也别一上来就装一堆用不上的东西环境越干净出问题时越容易定位。# 1) 确认 Python 与三个核心依赖齐全 python -c import numpy, scipy, matplotlib; print(deps ok) # 2) 进到解压后的项目目录看示例入口有哪些 ls examples/ # 3) 跑二维追逃那个示例文件名以解压后实际为准挑名字带 pursuit 的 python examples/pursuit_evasion_2d.py第一行只做导入检查不打印版本号因为版本兼容问题在数值包里通常表现为AttributeError而不是精度差异导入通过就够。第二行列目录是为了确认示例的组织方式有些包按博弈类型分文件有些按求解器分看清楚再跑省得在错误的入口上浪费半天。第三行是冒烟测试能出图就说明整条链路是通的后面再改自己的模型。如果第二步列出来的目录里只有求解器没有示例说明这个包是库形态需要自己写驱动脚本。这时候别急着从零写先把第 3.2 节的博弈定义类抄进去接上第 3.3 节的求解函数就能跑起来。3.2 定义一个追逃博弈动力学与代价函数我把博弈定义成一个类好处是动力学、代价、控制集绑在一起换问题只换类求解器代码一行不动。下面这个二维追逃定义包含捕获半径和单位时间惩罚是这类问题最常见的配置。import numpy as np class PursuitEvasion2D: 二维追逃追方最小化、逃方最大化代价为终端脱靶量 时间惩罚 state_dim 2 # 相对坐标 [rx, ry]已用平移不变性降维 u_dim 2 v_dim 2 def __init__(self, u_max1.0, v_max0.8, capture_radius0.1, time_penalty0.05): self.u_max u_max # 追方速度上限必须大于 v_max self.v_max v_max # 逃方速度上限 self.r_cap capture_radius # 捕获半径进入即视为捕获 self.c time_penalty # 每单位时间的代价越大越急于结束 def dynamics(self, t, r, u, v): # 相对位置导数追方靠近、逃方远离 return u - v def running_cost(self, t, r, u, v): # 常值时间惩罚不含控制能量项保持可追性判断清晰 return self.c def terminal_cost(self, r): # 捕获半径外按距离平方惩罚半径内为 0 d float(np.hypot(r[0], r[1])) return 1e3 * max(0.0, d - self.r_cap) ** 2dynamics返回的是相对位置速度之差这里的前提是双方都是简单积分器模型没有惯性、没有转向角约束。真实工程里如果换成有姿态动力学的模型这一行就要改成实际的状态方程同时状态维数会上升网格法可能就用不了了。time_penalty是权衡项。取值太小值函数在远离原点的地方几乎是平的梯度估计噪声会主导策略方向取值太大追方会表现为“急着结束”而在远处就开始猛冲轨迹看起来不自然。0.05 到 0.2 是个比较稳的区间具体按你的捕获半径和相对速度调。terminal_cost前面的 1e3 是量纲缩放目的是让终端项在总代价里占主导。如果你把它设成 1时间惩罚会在长轨迹上累积到和终端项同一量级值函数曲面会变得又缓又平梯度方向不可靠。这类权重我一般按“终端项典型值比积分项整段积分大一个量级”来定。3.3 在相对坐标网格上回推值函数有了博弈定义接下来在相对坐标平面上撒网格从终端时刻往回推进。核心更新式是 $V(t-\Delta t)V(t)\Delta t\cdot H$其中 $H(v_{max}-u_{max})|\nabla V|c$。下面是不带迎风修正的基础版本先跑通改进放到第 5 章讲。import numpy as np def solve_hji_relative(env, half2.0, nx161, ny161, T4.0, dt0.005): 在相对坐标网格上从终端时刻反推值函数 xs np.linspace(-half, half, nx) ys np.linspace(-half, half, ny) dx, dy xs[1] - xs[0], ys[1] - ys[0] X, Y np.meshgrid(xs, ys, indexingij) r0 np.stack([X, Y], axis0) # 终端条件捕获半径外的距离平方惩罚 V env.terminal_cost(r0) * np.ones_like(X) a env.v_max - env.u_max # 速度差负值表示追方占优 if a 0: raise ValueError(v_max u_max逃方永远逃得掉值函数无界) n_steps int(T / dt) for k in range(n_steps): Vy, Vx np.gradient(V, dy, dx) # 注意返回顺序是 [d/dy, d/dx] mag np.sqrt(Vx ** 2 Vy ** 2) V V dt * (a * mag env.c) # 时间反推的一阶推进 V np.maximum(V, 0.0) # 值函数非负抑制数值下冲 return xs, ys, Vnp.gradient(V, dy, dx)的返回顺序是[对第二维的导数, 对第一维的导数]因为网格是indexingij生成的第一维对应xs所以拿到的是Vy, Vx。这个顺序如果写反梯度方向会整体旋转 90 度策略看起来“有理有据”但完全错误属于典型的静默 bug。T的取值要覆盖你关心的最长决策时长。取得比实际捕获时间短值函数在远端的值还没传播到中心近端的策略是对的远端全是初值取得太长步数线性增加收益却很小。我的习惯是取理论捕获时间上界的 1.5 倍。nx、ny和half是一对需要平衡的量。half固定时格点数决定分辨率分辨率太粗梯度估计变成阶梯状策略会在格点之间来回横跳分辨率太细同样步数下推进成本上升而且边界外推误差会更明显。161 是个不错的起点先跑通再加密。3.4 从值函数梯度里抽出鞍点策略并仿真轨迹值函数本身只是一张表能执行的是策略。追方沿梯度下降方向走降低代价逃方沿梯度上升方向走提高代价绝对值都顶到各自的速度上限这就是鞍点策略的标准形式。import numpy as np def simulate_relative(env, xs, ys, V, r0(-1.6, 1.2), dt0.005, t_max12.0): dx, dy xs[1] - xs[0], ys[1] - ys[0] Vy, Vx np.gradient(V, dy, dx) def grad_at(p): # 最近邻查表生产环境建议换双线性插值 i int(np.clip(round((p[0] - xs[0]) / dx), 0, len(xs) - 1)) j int(np.clip(round((p[1] - ys[0]) / dy), 0, len(ys) - 1)) return Vx[i, j], Vy[i, j] r np.array(r0, dtypefloat) traj, t [r.copy()], 0.0 while t t_max and np.hypot(r[0], r[1]) env.r_cap: gx, gy grad_at(r) n np.hypot(gx, gy) 1e-12 # 防止除零远离原点时梯度基本不为零 u -env.u_max * np.array([gx, gy]) / n # 追方沿梯度反方向最小化 v env.v_max * np.array([gx, gy]) / n # 逃方沿梯度正方向最大化 r r dt * (u - v) traj.append(r.copy()) t dt return np.array(traj), t env PursuitEvasion2D(u_max1.0, v_max0.8) xs, ys, V solve_hji_relative(env, half2.0, nx161, ny161, T4.0, dt0.005) traj, t_cap simulate_relative(env, xs, ys, V) print(f捕获用时 {t_cap:.2f}s轨迹点数 {len(traj)}终点距原点 {np.hypot(*traj[-1]):.4f})梯度方向乘以负号还是正号取决于你把谁放在最小化位置。这里追方是最小化方所以沿 $-\nabla V$逃方是最大化方沿 $\nabla V$。这两个符号一旦反了你会看到一条“追方越追越远、逃方主动送上门”的诡异轨迹而代码本身不会报任何错。dt在这里是仿真步长和求解器里的 $dt$ 可以不同仿真步长只影响轨迹平滑度取小一点没关系。t_max是兜底防止策略有问题时死循环。最后打印的终点距原点应该略小于r_cap如果是 0.3 以上说明值函数还没推够步数或者网格分辨率不够。跑出来的轨迹形状通常是一条追方先斜插、逃方绕圈的曲线捕获时间随初始距离增长大致线性。如果轨迹是一条直线八成是逃方策略没生效如果轨迹在原点附近疯狂打转那是梯度在捕获半径附近趋零导致的属于正常现象加个半径内直接终止的判断就行。4. 求解器参数怎么设时间步长、网格分辨率与收敛判据同一份代码参数不一样结果能差出一倍。这一章讲三个必须自己定的量时间步长、空间分辨率、终止判据以及每个量的调节逻辑和观测指标。4.1 时间步长与显式推进的稳定性上限显式推进 HJI 方程是有硬稳定性上限的本质上是信息在网格上的传播速度不能超过一个格点每步。对当前这个模型传播速度上界是双方速度上限之和稳定性要求大致是$$\Delta t \le \frac{\min(\Delta x,\Delta y)}{\sqrt{2},(u_{max}v_{max})}$$按 $\Delta x0.025$、$u_{max}v_{max}1.8$ 算上限约 0.0098 秒。超过这个值值函数会出现高频振荡并逐步放大最终变成一片噪声。import numpy as np env PursuitEvasion2D(u_max1.0, v_max0.8) half, n 2.0, 161 dx 2 * half / (n - 1) dt_limit dx / (np.sqrt(2) * (env.u_max env.v_max)) print(f网格步长 {dx:.5f}稳定步长上限 {dt_limit:.5f}) for dt in [0.020, 0.010, 0.005, 0.002]: xs, ys, V solve_hji_relative(env, halfhalf, nxn, nyn, T4.0, dtdt) center V[n // 2, n // 2] rough np.abs(np.diff(V, axis0)).max() # 相邻格点最大跳变振荡时会暴涨 print(fdt{dt:.3f} 中心值{center:8.3f} 最大跳变{rough:8.3f})判断标准很简单dt超过上限时最大跳变会出现明显异常而中心值反而可能看起来正常——这就是为什么只看一个数值不够必须看局部粗糙度。我一般的习惯是取上限的三分之一到二分之一既不浪费步数也留出安全余量。4.2 空间分辨率与边界条件分辨率这块half和nx要一起调。判断够不够看两件事一是轨迹有没有贴到方框边界二是加密一倍后中心区域的值函数变化大不大。前者是覆盖问题后者是精度问题。边界条件是更容易被忽略的一项。上面代码里没有任何显式边界处理np.gradient在边界上用单侧差分等价于外推。对追逃问题正确做法通常在边界上把值函数固定成终端代价的初值或者干脆在边界附近加一个大惩罚让轨迹“不愿意”靠近边界。因为一旦轨迹真的贴边说明你的求解域取小了这时候边界处的结果本来就不可信。def solve_hji_with_boundary(env, half2.0, nx161, ny161, T4.0, dt0.005, boundary_penalty50.0): 边界格点强制回到初值避免轨迹贴边时用到不可信区域 xs np.linspace(-half, half, nx) ys np.linspace(-half, half, ny) dx, dy xs[1] - xs[0], ys[1] - ys[0] X, Y np.meshgrid(xs, ys, indexingij) V env.terminal_cost(np.stack([X, Y], axis0)) * np.ones_like(X) V0 V.copy() a env.v_max - env.u_max for _ in range(int(T / dt)): Vy, Vx np.gradient(V, dy, dx) V V dt * (a * np.sqrt(Vx ** 2 Vy ** 2) env.c) V np.maximum(V, 0.0) V[0, :] V0[0, :] boundary_penalty # 四条边固定并抬高 V[-1, :] V0[-1, :] boundary_penalty V[:, 0] V0[:, 0] boundary_penalty V[:, -1] V0[:, -1] boundary_penalty return xs, ys, Vboundary_penalty的作用是让边界在策略层面变成“高代价区”追方不愿去、逃方也不愿去从而把轨迹约束在内部可信区域。这个值不需要精确只要明显大于内部值函数的典型量级就行加太大会让梯度在边界附近出现陡壁插值误差变大。4.3 收敛阈值与终止判据时间推进法不像迭代法那样有天然的残差序列但可以自己构造两个监控量。第一个是相邻若干步之间值函数的最大变化第二个是值函数的变化率是否已经趋于常数。对追逃问题还有一个更强的判据跑一遍最优策略仿真看捕获时间是否随网格加密收敛。import numpy as np def solve_with_monitor(env, half2.0, nx161, ny161, T4.0, dt0.005, check_every200, tol1e-4): xs np.linspace(-half, half, nx) ys np.linspace(-half, half, ny) dx, dy xs[1] - xs[0], ys[1] - ys[0] X, Y np.meshgrid(xs, ys, indexingij) V env.terminal_cost(np.stack([X, Y], axis0)) * np.ones_like(X) a env.v_max - env.u_max n_steps int(T / dt) for k in range(n_steps): V_prev V Vy, Vx np.gradient(V, dy, dx) V V dt * (a * np.sqrt(Vx ** 2 Vy ** 2) env.c) V np.maximum(V, 0.0) if (k 1) % check_every 0: delta np.abs(V - V_prev).max() if delta tol: print(f第 {k1} 步后增量 {delta:.2e}低于阈值提前停止) break return xs, ys, Vcheck_every不要设成 1每步都算一次全网格差值在大的时候有明显开销。tol的物理意义是“再推一步值函数几乎不变”这时候再推就是在重复劳动。但要注意tol设得太松会在信息还没传播到中心区域时就停下表现为近原点处的策略正确、远处的策略乱指。稳妥做法是先不设阈值跑满 T记下每 200 步的增量曲线再根据曲线拐点定tol。4.4 参数速查表下表是我在二维追逃这类问题上用的默认值和调整方向不同问题形态会有偏移但相对关系基本成立。参数常用值调大的影响调小的影响该看什么指标dt稳定上限的 1/3出现高频振荡最大跳变暴涨步数线性增加收益递减相邻格点最大跳变nx、ny161梯度更平滑成本平方增长策略在格点间横跳加密一倍后中心值变化half典型捕获距离 × 2.5分辨率被稀释边界外推误差占比上升轨迹贴边界结果不可信轨迹最小边界间距T理论捕获时间 × 1.5步数增加远端收益很小中心区域策略未传播到位远端值函数是否仍是初值c0.050.2追方远处猛冲轨迹不自然值函数过平梯度噪声主导值函数曲面斜率u_max/v_max1.25 左右追方优势大值函数很快收敛逼近 1 时值函数无界风险中心值是否有限表里u_max/v_max那一行值得反复看。它不只是个数值本质是问题有没有解的前提。比值接近 1 时捕获时间对初始距离极其敏感值函数远端的值会变得非常大网格截断引入的相对误差被放大看起来像是求解失败其实是问题本身就在临界点上。5. 避坑与排查五条让值函数一夜之间变成噪声的记录数值求解这类方程出错方式非常隐蔽程序不崩、图也能画但结果就是不对。这一章按我实际踩过的顺序列五条每条都写清现象、原因和解决方式。5.1 值函数曲面出现锯齿甚至负值现象是画出来的值函数曲面像搓衣板沿着某个方向有一道道条纹捕获半径附近还冒出负值。第一反应通常是分辨率不够但其实加密到 401 之后锯齿只是变细了并没有消失。原因是中心差分格式用错了方向。HJI 方程属于双曲型信息沿特征方向单向传播用对称的中心差分会在梯度变化剧烈处产生非物理的振荡这在流体和光学计算里是老问题。np.gradient用的正是中心差分。解决办法是改成迎风差分根据系数符号选择前向或后向差分让信息只从上游取。当前模型系数 $av_{max}-u_{max}0$梯度模长应该用后向差分近似。改完之后锯齿基本消失代价是曲面在拐角处会略微变钝这是迎风格式的固有耗散。5.2 轨迹一路贴边滑到角落现象是仿真轨迹从初始点直接朝方框边缘走贴着边界滑到角上然后卡住打印出来的捕获时间等于t_max。原因有两个层面。表层原因是求解域half取小了实际捕获需要的横向距离超过了截断范围。深层原因是边界条件用了外推np.gradient在边缘单侧差分时等于默认“值函数继续线性延伸”策略会把轨迹往这个虚假的下降通道里推。解决分两步先把half放大到轨迹最小边界间距的 3 倍以上再给边界格点加固定惩罚值让边界成为高代价区。这两个动作要一起做只放大不处理边界轨迹不贴边了但边界附近的值函数仍然是外推出来的不可信。5.3 跑了十万步残差就是降不下去现象是监控日志里每 200 步的增量一直在同一量级徘徊中心值单调增长没有收敛迹象。先查两件事u_max是否大于v_max以及dt是否超过稳定上限。这两个都排除之后第三个可能是time_penalty太大。单位时间惩罚过大时值函数随时间近似线性增长增量永远降不下去因为代价本身在持续累积。解决办法是把c调小一个量级再跑一遍观察增量曲线是否出现拐点。如果调小后收敛了说明原来的c让积分项压过了终端项问题设定本身需要重新审视而不是求解器有问题。这个判断顺序很重要先排除设定错误再怀疑数值实现。5.4 策略在两个方向之间来回切换现象是轨迹出现密集的左右横跳看起来像抖动但值函数曲面本身很光滑。原因是梯度在某个区域内接近零除以模长时数值噪声被放大到满量程策略方向由噪声决定。常见于两个位置捕获半径附近的“谷底”以及值函数曲面斜率刚好为零的对称轴。解决办法是在策略提取时加一个方向死区模长小于阈值时保持上一步方向而不是重新计算。同时把最近邻查表换成双线性插值让梯度在格点之间连续变化抖动会明显减轻。这两个改动加起来不到十行代码但效果比调任何参数都直接。5.5 同样的脚本换台机器结果不一样现象是同一份代码、同一个随机种子如果有的话在另一台机器上跑出来的捕获时间差了好几个百分点。原因通常在浮点归约顺序上。如果代码里用了多线程的数组运算或者并行的 ODE 求解累加顺序不同会导致最低位差异经过几万步推进后被放大。另一种可能是不同机器的numpy底层 BLAS 实现不同。解决办法是排查阶段固定单线程、固定插值方式、去掉任何基于时间的种子先把可复现性拿到手。确认结果稳定之后再逐步放开性能选项每放开一个就对比一次数值看偏差是否在你关心的精度以内。这个习惯很啰嗦但在写进报告之前必须走一遍否则别人复现不出来你连解释的机会都没有。6. 进阶与验证怎么证明那张值函数图不是画得好看值函数曲面画出来是漂亮的漏斗形不说明它对。我在实际项目里被问过最多的问题就是“这个解可信吗”后来养成了一套三步验证的习惯每次出结论之前强制走一遍。第一步是回代残差。把求出来的策略代回 HJI 算子看每个网格点上的残差量级正常应该在 $10^{-3}$ 到 $10^{-2}$ 之间出现孤立的尖峰就说明那个区域格式有问题。第二步是策略扰动测试给逃方策略叠加一个小扰动统计平均捕获时间的变化理论上最优策略下捕获时间应该是所有扰动中的最小值如果随便扰动一下反而更快捕获说明你的策略根本没到鞍点。import numpy as np def perturb_test(env, xs, ys, V, r0(-1.6, 1.2), eps_list(0.0, 0.05, 0.1, 0.2), trials50): 给逃方方向叠加随机旋转扰动统计平均捕获时间 rng np.random.default_rng(7) dx, dy xs[1] - xs[0], ys[1] - ys[0] Vy, Vx np.gradient(V, dy, dx) def grad_at(p): i int(np.clip(round((p[0] - xs[0]) / dx), 0, len(xs) - 1)) j int(np.clip(round((p[1] - ys[0]) / dy), 0, len(ys) - 1)) return Vx[i, j], Vy[i, j] out {} for eps in eps_list: times [] for _ in range(trials): r np.array(r0, dtypefloat) t, dt 0.0, 0.005 while t 12.0 and np.hypot(r[0], r[1]) env.r_cap: gx, gy grad_at(r) n np.hypot(gx, gy) 1e-12 # 追方保持最优逃方方向加旋转扰动 theta rng.uniform(-eps, eps) gv np.array([gx, gy]) / n rot np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) u -env.u_max * gv v env.v_max * (rot gv) r r dt * (u - v) t dt times.append(t) out[eps] float(np.mean(times)) return out env PursuitEvasion2D(u_max1.0, v_max0.8) xs, ys, V solve_hji_relative(env, half2.0, nx161, ny161, T4.0, dt0.005) print(perturb_test(env, xs, ys, V))eps0.0那一档对应的就是最优逃方策略其它档位应该全部比它更大。如果某一档明显更小说明该区域的梯度方向不对多半是迎风方向或者边界处理的问题回到第 5 章逐条排查。trials取 50 是为了把随机扰动的影响平均掉太少会让结果抖得看不懂。第三步是网格收敛性。把nx从 81 加密到 161 再加密到 321记录同一个初始点下的捕获时间如果三次结果依次向某个值靠拢并且差值在成倍缩小说明离散化误差在正常收敛如果差值不缩反增那这张图再好看也不能用。这三步走完大概多花二十分钟但能省掉后面被人追着问“你怎么保证”的麻烦。从那以后我每次给对抗决策类的结果出图都强制先跑残差、扰动、加密这三遍宁可慢一点也不想在评审现场解释为什么曲线换台机器就变了。希望帮到你。本文还有配套的精品资源点击获取