PyBaMM 求解器收敛失败：先分清数值问题还是物理问题

一次长时间的 PyBaMM 退化模拟在第 300 圈崩掉，SUNDIALS 吐出一大段输出，末尾是 Corrector convergence failed repeatedly or with |h| = hmin 之类的字样。第一反应通常是去调容差或者把步长改小。这个反应基本上是错的，而且顺着它走会浪费掉几个小时。

这篇讲的是怎么区分两种情况。收敛失败要么是数值问题，要么是物理问题，而两者的修法方向相反。在一个物理上不可能的状态上放宽容差，救不回这次运行——它只是把”模型从哪里开始失去意义”这件事藏起来，然后让求解器吐出一批你之后会当真的数字。

一、这条报错到底在说什么

IDAKLU 是 BDF 求解器：每一步构造一个隐式代数方程组，用牛顿迭代去解。所谓 “corrector” 就是那个牛顿循环。有两件事会终结一步：

• corrector 不收敛。牛顿迭代跑到上限，残差仍然没降下去。求解器的反应是缩小步长重试。
• 步长撞到 hmin。反复缩小之后已经无可再缩。这就是 |h| = hmin 的含义。

关键推论就藏在这个顺序里：当你看到这条报错时，求解器已经把步长缩小过很多次了。手动调小步长或输出间隔，是在重复自适应控制器已经做过且失败了的事。如果缩步长有用，它早就有用了。

所以这条消息不是”步长太大”，它更接近于”在这个状态下，我愿意采用的任何步长都无法让牛顿迭代找到解”。

二、真正要先做的判断：数值还是物理

在动任何一个求解器选项之前，先确定你面对的是哪一种。这个测试很便宜。

import pybamm

# 用更简单的模型跑同一个实验。SPMe 放弃了 DFN 的空间分辨率但保留电解液梯度，
# 所以它能挺过一些 DFN 表达不了的状态。
for cls in (pybamm.lithium_ion.SPMe, pybamm.lithium_ion.DFN):
    sim = pybamm.Simulation(cls(), parameter_values=params, experiment=experiment)
    try:
        sol = sim.solve()
        print(f"{cls.__name__}: 完成 {len(sol.cycles)} 圈")
    except Exception as exc:
        print(f"{cls.__name__}: 失败 - {type(exc).__name__}: {exc}")

结果这样读：

实践中长时间退化模拟里第一行最常见，而它恰恰是人们用错工具、耗时最久的那一种。

三、找出那个物理上不可能的状态

当答案是物理问题，问题就变成了哪个变量越了界。嫌疑对象是那些方程默认会待在范围内的有界量：

• 颗粒表面浓度触到 0 或最大值。Butler-Volmer 动力学里含有 c_s_surf 和 (c_s_max - c_s_surf) 项，任一趋近于零，交换电流密度就塌掉、过电势就发散。
• 电解液浓度在隔膜或某个孔隙里降到零，电导率和扩散系数随之无定义。
• 孔隙率被 SEI 生长耗到零。这一条是长时间退化模拟特有的，也正是”模拟能跑完 250 圈却在第 251 圈失败”的原因。

看清楚它的方法是解到失败前再检查，而不是盯着求解器日志：

import numpy as np

# 先跑一段能成功的较短运行，再看有界变量离各自的极限有多近。
# 一次"跑通了但已经贴边"的运行，和几十圈之后的那次失败是同一回事。
sol = sim.solve()

for name in ("Negative particle surface concentration",
             "Positive particle surface concentration",
             "Electrolyte concentration [mol.m-3]",
             "Negative electrode porosity"):
    try:
        v = sol[name].entries
    except KeyError:
        continue          # 变量名随模型选项而变，缺了就跳过，别让整个诊断崩掉
    print(f"{name:<44} min={np.min(v):.4g}  max={np.max(v):.4g}")

化学计量数达到 0.999 或 0.001、电解液浓度比初值低一个数量级、孔隙率朝零走——这些都指向同一个结论：实验在要求电芯做它做不到的事。修法在协议或参数里，不在求解器里。

四、确实是数值问题时，有哪些选项

下面是 IDAKLU 文档记载的默认值。知道这些数字很重要，因为论坛上多数建议在调它们时并不说明原值是多少。

陷阱：设了却静默不生效的选项

IDAKLU 的 options 字典只在 model.convert_to_format == "casadi" 时生效。在其他格式的模型上设置它们，会被照单全收然后悄悄忽略。你会得出"这个选项没用"的结论，而它根本没被应用过。

solver = pybamm.IDAKLUSolver(
    rtol=1e-6,
    atol=1e-8,
    options={"max_nonlinear_iterations": 100, "max_order_bdf": 3},
)

model = pybamm.lithium_ion.DFN()

# 用断言而不是相信默认：这一行如果不是 "casadi"，上面那组 options 全是摆设
assert model.convert_to_format == "casadi", (
    f"convert_to_format 不是 'casadi' 时 IDAKLU 的 options 会被忽略；"
    f"当前是 {model.convert_to_format!r}"
)

sim = pybamm.Simulation(model, parameter_values=params,
                        experiment=experiment, solver=solver)

五、为什么长跑会失败而短跑不会

只在几百圈之后才出现的失败，几乎从来不是"第 3 圈还合适、第 300 圈就不合适"的求解器设置问题。一定有某种单调累积的东西。

机制通常是这样：退化子模型会改变方程所在的几何。SEI 生长消耗孔隙率，析锂和活性物质损失挪动可用的化学计量窗口。每一圈模型都离边界更近一点，最终一个原本只是刚性的状态变成了奇异的。

这对调试有两个直接推论：

• 用加速退化来更快复现，而不是跑更久。把 SEI 速率常数调高，直到同样的失败出现在第 20 圈，然后在那里诊断。二十圈的复现是你可以反复迭代的。
• 报告失败时一定要带上圈数。"参数集 X 下第 500 圈中的第 312 圈失败"是可诊断的，"模拟崩了"不是。

import pybamm

# 用加速退化把 300 圈才出现的失败压到 20 圈内复现 ——
# 目的不是得到有物理意义的结果，而是把迭代周期从几小时缩到几分钟
fast = params.copy()
fast["SEI kinetic rate constant [m.s-1]"] *= 50

sim = pybamm.Simulation(
    pybamm.lithium_ion.DFN(options={"SEI": "solvent-diffusion limited"}),
    parameter_values=fast,
    experiment=pybamm.Experiment(
        [("Discharge at 1C until 3.0V", "Charge at 1C until 4.2V", "Hold at 4.2V until C/50")] * 30
    ),
)
try:
    sol = sim.solve()
    print(f"跑完了 {len(sol.cycles)} 圈 —— 把倍数再调大")
except Exception as exc:
    print(f"复现成功: {type(exc).__name__}")

六、按这个顺序排查

1. 记录失败位置。圈数、圈内哪一步、最后一个成功的时间点。
2. 换 SPMe 重跑。如果同样失败，停止调求解器——问题在参数或协议。
3. 检查有界变量，在能成功的最长那次运行上看。找化学计量数贴近 0 或 1、电解液浓度塌陷、孔隙率趋零。
4. 把实验协议和电芯对一遍。退化后的电芯已经达不到的截止电压、已经撑不住的倍率，产生的正是这种失败。
5. 到这一步才轮到选项。先 max_nonlinear_iterations，再 max_order_bdf，而且在相信它生效之前先断言 convert_to_format == "casadi"。
6. 放宽过任何东西之后都要重新验证结果。靠放松接受阈值才跑完的运行，不等于可信的运行；拿一小段和更严格的设置对照一次。

七、跑通之后该记录什么

一次花掉一整天诊断的收敛失败，半年后会在另一组参数上重演。让第二次变便宜的是记录：

模型:     DFN, options {SEI: solvent-diffusion limited, ...}
参数:     OKane2022, 修改: SEI kinetic rate constant = 1e-15
实验:     1C/1C CCCV, 请求 500 圈
失败:     第 312 圈, 放电步, t = ...
诊断:     负极颗粒表面化学计量数 -> 0.9993
根因:     容量衰减约 30% 后截止电压已不可达
修法:     截止电压改到 2.8V；不是求解器改动
求解器:   IDAKLU, rtol 1e-4, atol 1e-6, 其余默认

注意"不是求解器改动"这一行。这份记录一半的价值在于拦住下一个人——通常就是你自己——再去动容差。

References

• PyBaMM IDAKLUSolver API 参考
• PyBaMM 讨论：corrector convergence failed repeatedly or with |h| = hmin
• PyBaMM 讨论：understanding solver doesn't converge
• PyBaMM 耦合退化 notebook

如果排查后发现是参数本身有问题而不是求解器，下一步见PyBaMM 参数拟合与可辨识性。

还可以参考：PyBaMM 架构与求解器选型（SPM/SPMe/DFN 的取舍，以及 IDAKLU 为什么适合长时间退化模拟）。