简介：difgames 是一个面向博弈论、控制理论与动态优化学习者的开源微分博弈项目，重点围绕哈密顿-雅可比-贝尔曼-伊萨克斯方程求解多决策者动态环境下的最优策略与纳什均衡，适合用于算法研究、课程实验和复杂决策问题探索。压缩包共 9 个文件，约 206KB，以 6 个 C 源码为主体，配套 1 个头文件、1 个 Shell 构建脚本和 1 份 PDF 文档，兼顾核心算法实现、图形显示、自动编译与理论说明。已有 240 人学习，读者可从中理解 HJBI 方程的数值近似思路，观察博弈过程的图形化结果与动态演化，并借助构建脚本在 Linux 环境快速编译运行。PDF 文档还可补充优化方法、算法原理与实验背景，为自动驾驶、经济模型或策略分析等场景提供可借鉴的开源实践样本，也适合教学演示、算法对照与二次开发。

1. 微分游戏的开源实现，为什么值得单独拆一遍

调一维追逃的参数调到凌晨两点，最磨人的往往不是数学推导，而是把耦合的两点边值问题塞进数值求解器之后那一串指数发散的残差。微分游戏（differential game）这套东西在教科书里推得很漂亮，落到代码上却常常卡在“初值怎么猜、Hamiltonian 怎么写、Riccati 边界条件给哪”这几个环节。difgames 这个开源包做的事情，就是把这些环节串成一条线：状态方程、代价函数、开环与反馈两类均衡、轨迹可视化，都放在同一套接口下。它适合做追逃对抗、多智能体博弈、制导律对比的研究生和算法工程师，也适合只想先跑通一个能收敛的例子、再回头补公式的人。下面按建模、求解、调参、排错、验证的顺序拆。

2. 微分博弈的数学骨架：从状态方程推到可运行的 Riccati

2.1 状态方程与代价函数：先把“谁在乎什么”写成矩阵

微分博弈的标准形式并不复杂：N 个玩家共享一个状态 x ∈ R^n，每人有自己的控制 u_i ∈ R^{m_i}，状态按 ẋ = f(x, u_1, …, u_N) 演化，第 i 个玩家最小化（或最大化）自己的代价

J_i = Φ_i(x(T)) + ∫₀ᵀ L_i(x, u_1, …, u_N) dt

纳什均衡的朴素定义是：任何人单独改自己的策略都占不到便宜。落到工程问题里，九成场景只需要两类——零和（追逃、拦截、制导对抗）和非零和（协同编队、资源分配）。零和的特殊之处在于 J_1 = -J_2，于是只需要解一个方程；非零和必须联立解两个。这个区别在代码里直接决定了 Riccati 是单组还是双组耦合，建模时就要定下来，别写到一半再改。

最容易翻车的地方是代价函数的物理含义。位置误差、控制能量、终端精度三者的权重比例，决定了游戏解出来是“追得凶”还是“省着动”。我一般会先把权重写成 dataclass，跑通再调，而不是散在若干函数里。

from dataclasses import dataclass, field import numpy as np @dataclass class LQGame: """线性二次微分博弈的最小描述，零和/非零和共用""" A: np.ndarray # 状态矩阵 (n, n) B: list # 每个玩家的控制矩阵 [(n, m_i), ...] Q: list # 状态权重，零和时两个玩家符号相反 R: list # 控制权重，必须正定 Qf: list = field(default=None) # 终端权重，None 表示不惩罚终端 T: float = 5.0 # 时域长度 x0: np.ndarray = field(default=None) def check(self): n = self.A.shape[0] assert self.A.shape == (n, n), "A 必须是方阵" for Bi, Qi, Ri in zip(self.B, self.Q, self.R): assert Bi.shape[0] == n, "B 的行数必须等于状态维数" assert np.allclose(Qi, Qi.T), "Q 必须对称" assert np.all(np.linalg.eigvalsh(Ri) > 1e-9), "R 必须正定" return self

这段check()看着啰嗦，但省掉的调试时间很值：R 不定时solve不会报错，只会悄悄给出一个符号错误的控制量，轨迹看上去“像那么回事”，结果完全反了。对称性检查同理，Q 写错一个非对角元，闭环特征值就可能跑到右半平面。

2.2 开环纳什与反馈纳什：两条路线怎么选

同样一个博弈，开环解和反馈解是两个不同的东西，很多初学者会把它们混为一谈。开环纳什把控制当成时间的函数 u_i(t)，所有人开局就把整条序列定死，数学上是一个两点边值问题（TPBVP）：状态方程正向、协态方程反向，边界条件一头给 x(0)、一头给协态 λ(T)。反馈纳什则把控制写成状态的函数 u_i = K_i(t)x，解的是耦合的 Hamilton-Jacobi-Bellman 方程组，在 LQ 情形退化成 Riccati 矩阵微分方程。

追逃这类问题我几乎只用反馈解。原因很简单：开环解是“我知道你会怎么跑，所以我这么追”，只要对方稍微偏一点，整条控制序列就废了；反馈解是“不管你跑到哪，我都按当前状态给加速度”，这才是工程上能用的东西。开环解仍然有价值——它可以作为反馈解的初值猜测，也可以用来校验反馈解在无扰动时是否退化到同一条轨迹。

选型上还有一条经验：状态维数 n ≤ 4、玩家数 2 的时候，直接上 Riccati 反向积分最快；n 上到两位数或者有控制约束，Riccati 就不够用了，得转向策略迭代或数值 HJB，这时候计算量会跳一个量级，别指望笔记本上几秒钟出结果。

2.3 零和 LQ 的 Riccati 方程：一次反向积分顶住大部分场景

零和 LQ 微分博弈值得单独写一遍，因为它是整个包里复用率最高的求解器。设定 J_1 被玩家 1 最小化、被玩家 2 最大化，值函数写成 V(x,t) = xᵀP(t)x，代入 HJB 方程并对 u_1 求最小、对 u_2 求最大，得到

Ṗ = -(Q + PA + AᵀP - P B₁ R₁⁻¹ B₁ᵀ P + P B₂ R₂⁻¹ B₂ᵀ P)，P(T) = Q_f

对应的最优反馈是 u_1 = -R₁⁻¹B₁ᵀPx、u_2 = +R₂⁻¹B₂ᵀPx。注意第二个式子的符号：玩家 2 是在最大化，梯度方向反过来，这个正负号写错是最高频的错误之一，现象是逃离者的控制量把距离越推越近。

import numpy as np from scipy.linalg import solve def riccati_rhs(t, p, A, B1, B2, Q, R1, R2): """零和 LQ 微分博弈的 Riccati 右端项，供 solve_ivp 反向积分""" n = A.shape[0] P = p.reshape(n, n) P = 0.5 * (P + P.T) # 每步强制对称，抑制长时程漂移 PB1 = P @ B1 PB2 = P @ B2 # 用 solve 而不是 inv：R 接近奇异时抛 LinAlgError，比静默算错好排查 u1_term = PB1 @ solve(R1, PB1.T) u2_term = PB2 @ solve(R2, PB2.T) dP = -(Q + P @ A + A.T @ P - u1_term + u2_term) return dP.ravel()

参数上有三处要留心。R1、R2是控制代价矩阵，数值越小代表该玩家越“舍得用力”，但太小会让solve接近奇异，一般保持在 1e-3 以上；Q在零和里两个玩家符号相反，写代码时只传一个正定的 Q，玩家 2 的符号由公式里的+ u2_term承担，不要重复取负；P(T) = Q_f是反向积分的起点，终端不惩罚就传零矩阵，别传 None。

非零和的情形需要两组耦合 Riccati 联立积分，未知量从 n² 变成 2n²，右端项里每个玩家的 P_i 都要用到对方的闭环矩阵。代码结构一样，只是把状态向量拼长，收敛性和稳定性会比零和差一些，时域长了容易积累误差，后面第 4 章会讲怎么压。

3. 跑通第一个例子：一维追逃从建模到出图

3.1 环境准备与依赖核对

这个包是纯 Python 实现，核心依赖只有 numpy 和 scipy，画图用 matplotlib，没有编译型依赖，装起来不折腾。我习惯用虚拟环境隔离，避免和系统里的 numpy 版本打架——scipy 对 numpy 的 ABI 比较敏感，版本错配的报错信息通常很难看懂。

python -m venv .venv source .venv/bin/activate # Windows 下用 .venv\Scripts\activate python -m pip install -U pip pip install "numpy>=1.24" "scipy>=1.11" "matplotlib>=3.7" # 从源码目录安装，走可编辑模式，改代码不用重装 cd difgames pip install -e .

装完先做一次导入自检，确认包路径指向的是源码目录而不是缓存里的旧版本：

python -c "import difgames, sys; print(difgames.__file__); print(sys.version)"

如果打印出来的路径在 site-packages 下而你确实在源码目录里改代码，说明装的是发布版，-e那步没生效，重来一次。这个坑我在两个不同项目上都踩过，症状是“明明改了代码，结果一点没变”，排查半小时才发现装的是另一份。

3.2 最小可运行示例：一维追逃的建模与求解

选一维追逃做第一个例子，是因为它状态只有两维、结果可以手算校验，跑通之后换成二维平面追逃只是把矩阵变大。状态取 x = [s, v]ᵀ，s 是追捕者与逃离者的相对位置，v 是相对速度；追捕者的加速度 u₁ 让 v 增大（追近），逃离者的加速度 u₂ 让 v 减小。

import numpy as np from scipy.integrate import solve_ivp # ---- 模型 ---- A = np.array([[0.0, 1.0], [0.0, 0.0]]) B1 = np.array([[0.0], [ 1.0]]) # 追捕者：正向加速度 B2 = np.array([[0.0], [-1.0]]) # 逃离者：反向加速度 Q = np.diag([1.0, 0.0]) # 只惩罚相对位置 R1 = np.array([[0.1]]) R2 = np.array([[0.1]]) Qf = np.diag([10.0, 0.0]) # 终端位置权重，逼它最后贴上去 T = 5.0 def riccati_rhs(t, p): n = 2 P = p.reshape(n, n) P = 0.5 * (P + P.T) PB1, PB2 = P @ B1, P @ B2 term = (Q + P @ A + A.T @ P - PB1 @ np.linalg.solve(R1, PB1.T) + PB2 @ np.linalg.solve(R2, PB2.T)) return (-term).ravel() # 反向积分：从 t=T 积到 t=0，初值是终端权重 Qf sol = solve_ivp(riccati_rhs, (T, 0.0), Qf.ravel(), method="Radau", rtol=1e-9, atol=1e-12, dense_output=True) assert sol.success, sol.message print("P(0) =", np.round(sol.sol(0.0).reshape(2, 2), 4))

这里method="Radau"不是随手选的。Riccati 方程反向积分时特征值可能很大，默认的 RK45 在时域稍长或权重稍极端时就会把步长压到极小，甚至报“步长过小”。Radau 是隐式方法，稳定性好，代价是每步要解线性系统，对这种 2×2 的小系统完全可以接受。rtol=1e-9偏严，是为了让 P 在接近终端时不出现明显的负特征值。

3.3 结果怎么读：轨迹、控制量与闭环稳定性

解出 P(t) 只是完成一半，真正要确认的是“这个博弈解是不是合理的”。我一般看三样东西：相对位置是否在终端附近收敛到零、两个玩家的控制量是否量级相当、闭环系统矩阵的特征值是否都在左半平面。

import matplotlib.pyplot as plt dt = 0.01 ts = np.arange(0.0, T + dt, dt) x = np.array([10.0, 0.0]) # 初始相对距离 10，相对速度 0 xs, u1s, u2s = [], [], [] for t in ts: P = sol.sol(t).reshape(2, 2) P = 0.5 * (P + P.T) u1 = -np.linalg.solve(R1, B1.T @ P) @ x u2 = np.linalg.solve(R2, B2.T @ P) @ x xs.append(x.copy()); u1s.append(u1[0]); u2s.append(u2[0]) x = x + dt * (A @ x + B1.ravel() * u1[0] + B2.ravel() * u2[0]) xs = np.array(xs) # 闭环特征值：任取中间时刻的 P 检查 Pm = sol.sol(T / 2).reshape(2, 2) K1 = -np.linalg.solve(R1, B1.T @ Pm) K2 = np.linalg.solve(R2, B2.T @ Pm) Acl = A + B1 @ K1 + B2 @ K2 print("闭环特征值:", np.linalg.eigvals(Acl)) print("终端相对位置:", xs[-1, 0])

读结果时有三条判据。第一，终端相对位置应该明显小于初值，但不必精确到零——终端权重 Q_f 是软约束，想要更贴近就把 Qf 加大，代价是控制量变大。第二，闭环特征值的实部必须为负，如果出现正实部，说明某个 R 或 Q 写反了符号，这时候轨迹往往表现为“越追越远”或者数值爆炸。第三，两个玩家的控制量曲线应该平滑，出现高频抖振通常不是博弈本身的性质，而是积分步长 dt 相对闭环带宽太大，把 dt 减半看抖振是否消失即可确认。

可视化部分把位置、速度、两组控制分别画三个子图，横轴是时间，用plt.subplots(3, 1, sharex=True)，比堆在一张图上好读得多。这一步没什么技术含量，但它是判断“解对不对”最直接的手段，别跳过。

4. 参数调优：时域、步长与收敛判据的实际取值

4.1 时域长度与积分步长怎么定

时域 T 是最容易被随手填的参数，填小了 P 还没收敛到稳态就结束，填大了浪费算力还放大数值误差。经验做法是先估一个闭环时间常数：把初始 P 代进去算一次瞬时增益，看闭环极点的时间常数 τ，取 T ≈ 4τ 到 6τ。追逃问题里 T 太短的典型现象是终端控制量突然翘起来——因为终端权重 Q_f 在最后一段主导了整条轨迹，博弈解退化成了“最后时刻才发力”。

这张表里的数值只解决“从哪开始”，具体问题还得扫一遍。我通常固定其他参数，把 q 按 0.1、1、10 三档扫，看终端误差和控制峰值的变化曲线，拐点附近就是可用的工作点。别用最优化算法自动搜——参数少的时候人工扫比写搜索脚本快，而且能看清每个参数在物理上干了什么。

4.2 收敛判据：看残差，不看“跑完了没”

solve_ivp返回success=True只代表积分器没报错，不代表解可用。Riccati 反向积分的有效性要额外验证两件事：P(t) 在整条时域上保持对称半正定，以及闭环系统在整个时域上稳定。前者可以用采样点上的最小特征值检查，后者在每个采样点算一次 Acl 的特征值实部最大值。

bad = [] for t in np.linspace(0.0, T, 101): P = sol.sol(t).reshape(2, 2) P = 0.5 * (P + P.T) if np.linalg.eigvalsh(P).min() < -1e-8: bad.append(("非半正定", t)) K1 = -np.linalg.solve(R1, B1.T @ P) K2 = np.linalg.solve(R2, B2.T @ P) if np.linalg.eigvals(A + B1 @ K1 + B2 @ K2).max().real > 1e-6: bad.append(("闭环不稳定", t)) print("异常采样点:", bad[:5], "共", len(bad))

这段检查比我见过的任何“收敛曲线”都管用。Riccati 方程本身对初值不敏感（初值就是 Q_f，确定的），但对数值误差敏感，误差累积到一定程度会先表现为 P 的最小特征值轻微变负，再演化成轨迹发散。早发现这一步，能省掉后面大量“为什么轨迹不对”的排查。

4.3 数值稳定性：从显式积分换到隐式，从单打靶换到松弛

稳定性问题分两层。第一层是积分器层面：Riccati 反向积分等价于在时间上倒着解一个可能不稳定的系统，正向稳定的特征值反过来就是快速增长的，显式方法必须把步长压到极小才能稳。换 Radau 或 BDF 是最省事的解法，代价是每步的线性代数开销，但对 n ≤ 10 的矩阵完全无所谓。

第二层是问题层面：如果时域很长、权重差异很大（比如 q=1 而 r=1e-4），Riccati 的解会在某些区间变化极快，任何固定容差都难受。这时候两个办法：一是把时域切成若干段，每段单独积分，用上一段的终点作为下一段的起点，段内用更细的容差；二是干脆放弃直接积分，改用策略迭代——先猜一组增益 K，解 Lyapunov 方程得到代价，再更新 K，反复直到收敛。策略迭代对初值不敏感，缺点是每轮要解一个 Lyapunov 方程，n 大了就慢。

我一般先用直接积分，出现“P 在中段剧烈震荡”再切策略迭代。切换前先确认不是 R 接近奇异导致的——把 R 放大 10 倍看震荡是否消失，是最快的判别手段。

5. 避坑清单：五个把我卡到下班的现场记录

5.1 环境与入口相关的两个坑

现象：在源码目录下执行python demo.py，报ModuleNotFoundError: No module named 'difgames'，但明明包就在旁边。

原因：仓库用的是 src 布局，包实际在src/difgames/下，Python 的模块搜索路径里没有src。直接跑脚本时当前目录被加进sys.path，但src不在其中。

解决：优先用pip install -e .走可编辑安装，把包注册进环境；临时调试可以PYTHONPATH=src python demo.py。别用sys.path.append硬改，换台机器就失效。

现象：改了源码里的求解器，重跑结果一点变化没有。

原因：之前装过一次发布版，import时优先命中了 site-packages 里的那份，源码目录只是“在旁边看着”。

解决：python -c "import difgames; print(difgames.__file__)"确认路径。如果指向 site-packages，先pip uninstall difgames再pip install -e .。这个坑很隐蔽，因为代码没报错，只是“静默地用错了版本”。

5.2 数值与结果解释相关的三个坑

现象：开环打靶法跑几十次，初值稍改一点就发散，残差从 1e-2 直接跳到 1e3。

原因：开环 TPBVP 的协态方程反向积分时存在快速增长的模态，单打靶对这个模态极其敏感，属于方法本身的病态，不是代码写错了。

解决：换多重打靶（把时域分 5~10 段，段间连续性作为额外约束），或者用反馈 Riccati 解生成一条标称轨迹，再拿它当打靶的初值。后者在实践中更省事，也顺便验证了两条求解路线的一致性。

现象：Riccati 积分到后半段，P 出现 NaN，solve_ivp返回success=False。

原因：多半是某个 R 接近奇异，或者容差太松导致误差在刚度大的区间被放大。也有可能是 Q 不对称，导致 P 的对称性被破坏后负特征值滚雪球。

解决：先np.linalg.eigvalsh(R)检查最小特征值，低于 1e-3 就把 R 放大；再确认 Q 严格对称（np.allclose(Q, Q.T)）；最后把rtol收紧到 1e-10、换 Radau。三步走完还 NaN，基本可以断定是模型本身在这个时域上无解，缩短 T 试一次。

现象：逃离者的控制量曲线越跑越大，相对距离不降反升，但代码没报任何错。

原因：Riccati 公式里玩家 2 的符号写反了。零和博弈中玩家 2 是最大化代价，控制律是+R₂⁻¹B₂ᵀPx，写成负号等于让它配合追捕者。

解决：用第 3 章那段闭环特征值检查，符号错的时候Acl大概率有正实部。另一个快速判别法是设 u₂ ≡ 0 跑一遍，看距离是否单调下降——如果关闭逃离者后反而收敛得更慢，符号基本可以确定有问题。

现象：策略迭代跑起来代价震荡，两个玩家的增益此消彼长，几十轮都不收敛。

原因：同步更新——两个玩家都用上一轮的对方策略算自己的最优响应，等价于在鞍点上做对称迭代，容易在两个状态之间来回跳。

解决：改成异步（Gauss-Seidel）更新，先更新玩家 1 再用新的玩家 1 更新玩家 2；或者给更新量加一个阻尼系数（新增益 = 0.5 旧 + 0.5 新）。阻尼系数从 0.5 起调，收敛慢但稳。

6. 进阶：换掉内置算例，用三重校验确认解是对的

内置算例跑通只说明环境和接口没问题，真正要用到自己的问题上，第一步是把动力学和代价换掉。以二维平面追逃为例，状态从 2 维扩到 4 维（两个方向的相对位置和相对速度），A 变成分块对角，B₁、B₂ 各扩成 4×2，Q 和 R 相应升维。代码结构一行不用改，只需要重新检查 R 的正定性和 Q 的对称性——维数一高，手写矩阵最容易在这两处出错。

换完之后别急着看轨迹好不好看，先做三重校验：

第三项最容易被忽略，但它能抓出绝大部分符号和权重错误。把 u₂ 强制置零，博弈退化成单人 LQR，此时用scipy.linalg.solve_continuous_are解代数 Riccati 方程可以得到解析最优增益，两条轨迹应当完全重合。我有一回就是因为这一步没做，把 Q 的符号写反了，跑出来的轨迹“看着挺合理”，直到写报告对比基线时才发现追捕者一直在帮逃离者加速。

还有一个小技巧：把时域按 5 段切开，每段单独积分，段间用上一段的 P 作为下一段的初值，最后拼起来和一次性积分的结果对比。差异超过 1e-3 就说明刚度太大，一次性积分不可信，得改用分段或策略迭代。这个习惯是从一次长时域制导仿真里逼出来的，那次 P 在最后 20% 的时域上偷偷变负，轨迹末端偏了将近一倍，肉眼完全看不出来。

从那以后我每次换新模型，都强制走一遍“关掉一个玩家对比 LQR”这一步，不管多赶时间。希望帮到你。

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