企业官网建设流程全解析

简介：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 f+L}=0$，终端条件 $V(T,x)=\varphi(x)$。这一串符号落到工程上就是四个必须显式定义的对象，缺一个都跑不起来。

控制集那一行是最容易翻车的地方。追方速度上限只要不大于逃方，在无界平面上就是永远追不上，代码不会报错，它只会让值函数在几万步迭代里线性地往上爬，你以为是在收敛，其实是在发散。我在第一次跑的时候就吃了这个亏，日志里残差一直降不下去，后来把 $u_{max}$ 从 0.8 调到 1.0，三分钟就收敛了。

2.2 相对坐标降维：四维状态怎么压成两维网格

追逃问题的自然状态是四个量：追方位置两维、逃方位置两维。四维网格求解的代价是灾难性的，161 个点每维就是 $161^4\approx6.7$ 亿个格点，普通机器连一个时间步都扛不住。但这类问题几乎都有平移不变性，把状态换成相对位置 $r=p_{追}-p_{逃}$，动力学变成 $\dot{r}=u-v$，状态只剩两维，运行代价直接掉到 $161^2\approx2.6$ 万个格点，单机几秒钟一步。

这个变换不是免费的。它要求代价函数只依赖相对位置和相对速度，一旦你的问题里有绝对位置相关的项（比如边界约束、地形代价、某个固定目标点），降维就不成立，得回到全状态空间或者用别的分解方式。判断方法很朴素：把动力学和代价里的绝对坐标全部替换成相对坐标，如果表达式中不再出现单独的 $p_{追}$ 或 $p_{逃}$，就可以降。

降维之后还有一个隐含假设：相对坐标平面是无限的。实际实现必须截断成一个有界方框，然后在边界上做处理。这个截断大小直接决定你的求解域够不够覆盖实际轨迹，取小了轨迹会贴边走，取大了网格分辨率被稀释，梯度估计噪声变大。我一般的做法是先用解析上界估一个粗的捕获距离，再乘 2.5 倍当半边长。

2.3 三条求解路线怎么选：矩阵方程、网格 PDE、迭代最佳响应

微分博弈的数值解法不是只有一条路，选错路线比调错参数更浪费时间。下面这张表是我自己在不同项目里踩过之后总结的选型依据，覆盖了绝大多数中小规模场景。

线性二次型那条路值得单独提一句。它把 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_max=1.0, v_max=0.8, capture_radius=0.1, time_penalty=0.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) ** 2

dynamics返回的是相对位置速度之差，这里的前提是双方都是简单积分器模型，没有惯性、没有转向角约束。真实工程里如果换成有姿态动力学的模型，这一行就要改成实际的状态方程，同时状态维数会上升，网格法可能就用不了了。

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, half=2.0, nx=161, ny=161, T=4.0, dt=0.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, indexing="ij") r0 = np.stack([X, Y], axis=0) # 终端条件：捕获半径外的距离平方惩罚 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, V

np.gradient(V, dy, dx)的返回顺序是[对第二维的导数, 对第一维的导数]，因为网格是indexing="ij"生成的，第一维对应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), dt=0.005, t_max=12.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, dtype=float) 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_max=1.0, v_max=0.8) xs, ys, V = solve_hji_relative(env, half=2.0, nx=161, ny=161, T=4.0, dt=0.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 x=0.025$、$u_{max}+v_{max}=1.8$ 算，上限约 0.0098 秒。超过这个值，值函数会出现高频振荡并逐步放大，最终变成一片噪声。

import numpy as np env = PursuitEvasion2D(u_max=1.0, v_max=0.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, half=half, nx=n, ny=n, T=4.0, dt=dt) center = V[n // 2, n // 2] rough = np.abs(np.diff(V, axis=0)).max() # 相邻格点最大跳变，振荡时会暴涨 print(f"dt={dt:.3f} 中心值={center:8.3f} 最大跳变={rough:8.3f}")

判断标准很简单：dt超过上限时，最大跳变会出现明显异常，而中心值反而可能看起来正常——这就是为什么只看一个数值不够，必须看局部粗糙度。我一般的习惯是取上限的三分之一到二分之一，既不浪费步数也留出安全余量。

4.2 空间分辨率与边界条件

分辨率这块，half和nx要一起调。判断够不够，看两件事：一是轨迹有没有贴到方框边界，二是加密一倍后中心区域的值函数变化大不大。前者是覆盖问题，后者是精度问题。

边界条件是更容易被忽略的一项。上面代码里没有任何显式边界处理，np.gradient在边界上用单侧差分，等价于外推。对追逃问题，正确做法通常在边界上把值函数固定成终端代价的初值，或者干脆在边界附近加一个大惩罚，让轨迹“不愿意”靠近边界。因为一旦轨迹真的贴边，说明你的求解域取小了，这时候边界处的结果本来就不可信。

def solve_hji_with_boundary(env, half=2.0, nx=161, ny=161, T=4.0, dt=0.005, boundary_penalty=50.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, indexing="ij") V = env.terminal_cost(np.stack([X, Y], axis=0)) * 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, V

boundary_penalty的作用是让边界在策略层面变成“高代价区”，追方不愿去、逃方也不愿去，从而把轨迹约束在内部可信区域。这个值不需要精确，只要明显大于内部值函数的典型量级就行，加太大会让梯度在边界附近出现陡壁，插值误差变大。

4.3 收敛阈值与终止判据

时间推进法不像迭代法那样有天然的残差序列，但可以自己构造两个监控量。第一个是相邻若干步之间值函数的最大变化，第二个是值函数的变化率是否已经趋于常数。对追逃问题，还有一个更强的判据：跑一遍最优策略仿真，看捕获时间是否随网格加密收敛。

import numpy as np def solve_with_monitor(env, half=2.0, nx=161, ny=161, T=4.0, dt=0.005, check_every=200, tol=1e-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, indexing="ij") V = env.terminal_cost(np.stack([X, Y], axis=0)) * 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"第 {k+1} 步后增量 {delta:.2e}，低于阈值，提前停止") break return xs, ys, V

check_every不要设成 1，每步都算一次全网格差值在大的时候有明显开销。tol的物理意义是“再推一步值函数几乎不变”，这时候再推就是在重复劳动。但要注意，tol设得太松会在信息还没传播到中心区域时就停下，表现为近原点处的策略正确、远处的策略乱指。稳妥做法是先不设阈值跑满 T，记下每 200 步的增量曲线，再根据曲线拐点定tol。

4.4 参数速查表

下表是我在二维追逃这类问题上用的默认值和调整方向，不同问题形态会有偏移，但相对关系基本成立。

表里u_max/v_max那一行值得反复看。它不只是个数值，本质是问题有没有解的前提。比值接近 1 时，捕获时间对初始距离极其敏感，值函数远端的值会变得非常大，网格截断引入的相对误差被放大，看起来像是求解失败，其实是问题本身就在临界点上。

5. 避坑与排查：五条让值函数一夜之间变成噪声的记录

数值求解这类方程，出错方式非常隐蔽：程序不崩、图也能画，但结果就是不对。这一章按我实际踩过的顺序列五条，每条都写清现象、原因和解决方式。

5.1 值函数曲面出现锯齿甚至负值

现象是画出来的值函数曲面像搓衣板，沿着某个方向有一道道条纹，捕获半径附近还冒出负值。第一反应通常是分辨率不够，但其实加密到 401 之后锯齿只是变细了，并没有消失。

原因是中心差分格式用错了方向。HJI 方程属于双曲型，信息沿特征方向单向传播，用对称的中心差分会在梯度变化剧烈处产生非物理的振荡，这在流体和光学计算里是老问题。np.gradient用的正是中心差分。

解决办法是改成迎风差分：根据系数符号选择前向或后向差分，让信息只从上游取。当前模型系数 $a=v_{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), trials=50): """给逃方方向叠加随机旋转扰动，统计平均捕获时间""" 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, dtype=float) 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_max=1.0, v_max=0.8) xs, ys, V = solve_hji_relative(env, half=2.0, nx=161, ny=161, T=4.0, dt=0.005) print(perturb_test(env, xs, ys, V))

eps=0.0那一档对应的就是最优逃方策略，其它档位应该全部比它更大。如果某一档明显更小，说明该区域的梯度方向不对，多半是迎风方向或者边界处理的问题，回到第 5 章逐条排查。trials取 50 是为了把随机扰动的影响平均掉，太少会让结果抖得看不懂。

第三步是网格收敛性。把nx从 81 加密到 161 再加密到 321，记录同一个初始点下的捕获时间，如果三次结果依次向某个值靠拢，并且差值在成倍缩小，说明离散化误差在正常收敛；如果差值不缩反增，那这张图再好看也不能用。

这三步走完大概多花二十分钟，但能省掉后面被人追着问“你怎么保证”的麻烦。从那以后我每次给对抗决策类的结果出图，都强制先跑残差、扰动、加密这三遍，宁可慢一点也不想在评审现场解释为什么曲线换台机器就变了。希望帮到你。

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