简介difgames 是一个面向博弈论、控制理论与动态优化学习者的开源微分博弈项目重点围绕哈密顿-雅可比-贝尔曼-伊萨克斯方程求解多决策者动态环境下的最优策略与纳什均衡适合用于算法研究、课程实验和复杂决策问题探索。压缩包共 9 个文件约 206KB以 6 个 C 源码为主体配套 1 个头文件、1 个 Shell 构建脚本和 1 份 PDF 文档兼顾核心算法实现、图形显示、自动编译与理论说明。已有 240 人学习读者可从中理解 HJBI 方程的数值近似思路观察博弈过程的图形化结果与动态演化并借助构建脚本在 Linux 环境快速编译运行。PDF 文档还可补充优化方法、算法原理与实验背景为自动驾驶、经济模型或策略分析等场景提供可借鉴的开源实践样本也适合教学演示、算法对照与二次开发。1. 微分游戏的开源实现为什么值得单独拆一遍调一维追逃的参数调到凌晨两点最磨人的往往不是数学推导而是把耦合的两点边值问题塞进数值求解器之后那一串指数发散的残差。微分游戏differential game这套东西在教科书里推得很漂亮落到代码上却常常卡在“初值怎么猜、Hamiltonian 怎么写、Riccati 边界条件给哪”这几个环节。difgames 这个开源包做的事情就是把这些环节串成一条线状态方程、代价函数、开环与反馈两类均衡、轨迹可视化都放在同一套接口下。它适合做追逃对抗、多智能体博弈、制导律对比的研究生和算法工程师也适合只想先跑通一个能收敛的例子、再回头补公式的人。下面按建模、求解、调参、排错、验证的顺序拆。2. 微分博弈的数学骨架从状态方程推到可运行的 Riccati2.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(defaultNone) # 终端权重None 表示不惩罚终端 T: float 5.0 # 时域长度 x0: np.ndarray field(defaultNone) 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 矩阵微分方程。维度开环纳什反馈纳什决策变量u_i(t)与状态无关K_i(t)随状态在线变化数学形式两点边值问题 TPBVP耦合 HJB / Riccati数值方法单打靶、多重打靶、松弛法、配点法Riccati 反向积分、策略迭代抗扰动差偏离标称轨迹后策略失效好状态偏离后控制自动修正计算量低未知量是控制序列中未知量是 n×n 矩阵典型场景离线规划、开环轨迹校验在线反馈、制导律、追逃对抗追逃这类问题我几乎只用反馈解。原因很简单开环解是“我知道你会怎么跑所以我这么追”只要对方稍微偏一点整条控制序列就废了反馈解是“不管你跑到哪我都按当前状态给加速度”这才是工程上能用的东西。开环解仍然有价值——它可以作为反馈解的初值猜测也可以用来校验反馈解在无扰动时是否退化到同一条轨迹。选型上还有一条经验状态维数 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 而不是 invR 接近奇异时抛 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 numpy1.24 scipy1.11 matplotlib3.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() # 反向积分从 tT 积到 t0初值是终端权重 Qf sol solve_ivp(riccati_rhs, (T, 0.0), Qf.ravel(), methodRadau, rtol1e-9, atol1e-12, dense_outputTrue) assert sol.success, sol.message print(P(0) , np.round(sol.sol(0.0).reshape(2, 2), 4))这里methodRadau不是随手选的。Riccati 方程反向积分时特征值可能很大默认的 RK45 在时域稍长或权重稍极端时就会把步长压到极小甚至报“步长过小”。Radau 是隐式方法稳定性好代价是每步要解线性系统对这种 2×2 的小系统完全可以接受。rtol1e-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, sharexTrue)比堆在一张图上好读得多。这一步没什么技术含量但它是判断“解对不对”最直接的手段别跳过。4. 参数调优时域、步长与收敛判据的实际取值4.1 时域长度与积分步长怎么定时域 T 是最容易被随手填的参数填小了 P 还没收敛到稳态就结束填大了浪费算力还放大数值误差。经验做法是先估一个闭环时间常数把初始 P 代进去算一次瞬时增益看闭环极点的时间常数 τ取 T ≈ 4τ 到 6τ。追逃问题里 T 太短的典型现象是终端控制量突然翘起来——因为终端权重 Q_f 在最后一段主导了整条轨迹博弈解退化成了“最后时刻才发力”。参数建议起点作用方向调错的典型现象时域 T4~6 倍闭环时间常数太短终端主导太长误差累积控制量末端翘起或前段几乎为零积分容差 rtol1e-8 ~ 1e-9越严越稳、越慢过松时 P 出现负特征值、轨迹后期发散正向步长 dt1e-2 ~ 1e-3需小于最快闭环带宽的 1/10控制量高频抖振、欧拉积分发散位置权重 q先取 1与 r 成 10:1 量级越大追得越紧过大导致控制饱和、数值刚度上升控制权重 r先取 0.1别低于 1e-3越小越舍得用力过小使solve接近奇异、控制量爆掉终端权重 Q_f初始状态权重的 5~10 倍越大末端贴得越紧过大导致末端控制量尖峰这张表里的数值只解决“从哪开始”具体问题还得扫一遍。我通常固定其他参数把 q 按 0.1、1、10 三档扫看终端误差和控制峰值的变化曲线拐点附近就是可用的工作点。别用最优化算法自动搜——参数少的时候人工扫比写搜索脚本快而且能看清每个参数在物理上干了什么。4.2 收敛判据看残差不看“跑完了没”solve_ivp返回successTrue只代表积分器没报错不代表解可用。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 的矩阵完全无所谓。第二层是问题层面如果时域很长、权重差异很大比如 q1 而 r1e-4Riccati 的解会在某些区间变化极快任何固定容差都难受。这时候两个办法一是把时域切成若干段每段单独积分用上一段的终点作为下一段的起点段内用更细的容差二是干脆放弃直接积分改用策略迭代——先猜一组增益 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 .走可编辑安装把包注册进环境临时调试可以PYTHONPATHsrc 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 出现 NaNsolve_ivp返回successFalse。原因多半是某个 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×2Q 和 R 相应升维。代码结构一行不用改只需要重新检查 R 的正定性和 Q 的对称性——维数一高手写矩阵最容易在这两处出错。换完之后别急着看轨迹好不好看先做三重校验校验项做法通过标准时间步长收敛dt 减半重跑比较终端代价两次差异 1%容差收敛rtol 收紧一个量级重跑P(0) 的元素变化 1e-4退化一致性关掉玩家 2 的控制与单方最优控制对比轨迹重合到 1e-3第三项最容易被忽略但它能抓出绝大部分符号和权重错误。把 u₂ 强制置零博弈退化成单人 LQR此时用scipy.linalg.solve_continuous_are解代数 Riccati 方程可以得到解析最优增益两条轨迹应当完全重合。我有一回就是因为这一步没做把 Q 的符号写反了跑出来的轨迹“看着挺合理”直到写报告对比基线时才发现追捕者一直在帮逃离者加速。还有一个小技巧把时域按 5 段切开每段单独积分段间用上一段的 P 作为下一段的初值最后拼起来和一次性积分的结果对比。差异超过 1e-3 就说明刚度太大一次性积分不可信得改用分段或策略迭代。这个习惯是从一次长时域制导仿真里逼出来的那次 P 在最后 20% 的时域上偷偷变负轨迹末端偏了将近一倍肉眼完全看不出来。从那以后我每次换新模型都强制走一遍“关掉一个玩家对比 LQR”这一步不管多赶时间。希望帮到你。本文还有配套的精品资源点击获取
)
