Most tutorials that treat PyBaMM as a tool start by running a discharge curve. This one goes the other way: it looks at how PyBaMM compiles electrochemical PDEs into a solvable system in the first place. The reason is practical — when the solver dies at cycle 300 in the middle of the night, you need to know which layer it died in, and that only comes from understanding the architecture.
One spelling note for searchers: the formal name is PyBaMM (Python Battery Mathematical Modelling), not PyBAMM. EIS used to live in a separate package, pybammeis, and has since merged into core — use pybamm.EISSimulation, which linearises the underlying DAE system in the frequency domain.
PyBaMM battery modelling series (4 parts): ① Architecture and solvers → ② EIS labels → ③ Dataset pipeline → ④ Training SOH and RUL. This is part 1.
1. PyBaMM as a DAE compilation pipeline
The critical architectural decision in PyBaMM is that electrochemistry is abstracted into a symbolic expression tree before any numerical discretisation. The Simulation class acts as a compiler front end, translating continuum equations into a discretised state space.
The instantiation order is not arbitrary — it is the compilation stages:
- Model: SPM, SPMe, or DFN built on P2D porous-electrode theory.
- Submodel options: SEI growth, plating overpotential, particle fracture mechanics, LAM, surface-area formulation. This must happen at model construction — pushing values into parameters afterwards enables nothing.
- Parameterisation: inject nonlinear functions such as OCP curves plus scalar properties, mapping symbols to physical quantities.
- Experiment: current density input, cutoff voltages, CCCV cycling, boundary conditions.
- Discretisation and solver: finite volume in space, PDE to DAE, handed to a BDF solver that can survive extreme stiffness.
Knowing this order has direct debugging value: the stage an error surfaces in tells you where to look. Missing parameters fail at stage 3, mesh problems at stage 5 — and “the model ran but nothing degraded” raises nothing at all, because it is a silent failure at stage 2.
2. Formalising DFN, and where stiffness comes from
DFN solves coupled conservation laws in the solid and electrolyte phases. Solid lithium concentration $c_s(r,x,t)$ inside a particle obeys Fick’s second law in spherical coordinates:
$$ \frac{\partial c_s}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 D_s(c_s) \frac{\partial c_s}{\partial r} \right) $$
This couples to the particle-electrolyte interface through Butler-Volmer kinetics, which sets the volumetric transfer current density $j(x,t)$:
$$ j = a_s i_0 \left[ \exp\left(\frac{\alpha_a F \eta}{R T}\right) – \exp\left(-\frac{\alpha_c F \eta}{R T}\right) \right] $$
Exchange current density $i_0$ depends on $c_s$, electrolyte concentration $c_e$ and local overpotential $\eta = \phi_s – \phi_e – U_{OCP}(c_s)$.
Stiffness can be pinned down more precisely than “the system is stiff”: characteristic time constants span roughly ten orders of magnitude — double-layer capacitance responds in microseconds, solid diffusion over hours. The stiffness ratio is that span. Explicit methods such as RK4 are bound by the CFL condition to steps short enough to resolve the fastest process, so simulating one hour would take on the order of a hundred million steps. That is not slow; it is infeasible.
3. Symbolic trees and automatic differentiation
PyBaMM does not hard-code sparse matrices. It builds equations from symbolic nodes (pybamm.Variable, pybamm.grad, pybamm.div). The snippet below runs as written and shows the expression tree and symbolic differentiation:
import pybamm
x = pybamm.Variable("x")
expr = 3 * x**2 + 2 * x
print(expr) # the expression tree itself
print(expr.diff(x)) # symbolic derivative - a new tree, not a number
A real diffusion operator looks like this. Note that grad and div require a variable with a domain — spatial operators need to know which geometry they discretise over, and omitting it raises immediately:
c_s = pybamm.Variable("Solid concentration", domain="negative particle")
D_s = pybamm.Parameter("Solid diffusion coefficient")
N_s = -D_s * pybamm.grad(c_s) # diffusive flux
dcdt = -pybamm.div(N_s) # rate of change of concentration
The symbolic graph is what lets PyBaMM swap parameter sets, discretise on any mesh, and use CasADi for automatic differentiation to build an exact analytical Jacobian. An exact Jacobian directly reduces Newton iterations per implicit step — not a nicety, but the difference between converging and not on a stiff system.
4. Choosing a solver
After spatial discretisation, DFN yields a large stiff system of ODEs plus algebraic constraints. PyBaMM compiles the symbolic tree into CasADi SX/MX graphs, generating optimised C code to evaluate the right-hand side and Jacobian, then hands the problem to SUNDIALS (IDA/IDAS or CVODES) and its variable-order, variable-step BDF methods. Step size adapts: tiny through a current step, minutes long during rest.
| Solver | Use for | Caveat |
|---|---|---|
CasadiSolver |
Default, most general | Noticeably slower than IDAKLU on long cycling experiments |
IDAKLUSolver |
First choice for long degradation runs | Sparse linear algebra; the advantage grows with problem size |
ScipySolver |
Quick experiments on pure-ODE models | Cannot handle DAEs with algebraic constraints |
Tolerance is another knob that must be managed explicitly. Defaults are fine on healthy cells, but in deeply degraded states they admit numerical noise that is indistinguishable from electrochemistry — you think you are seeing a degradation signature when the solver is simply jittering. For degradation work, record tolerance and solver version in the sample metadata.
5. The three models are different state spaces, not accuracy tiers
For machine learning practitioners this is the key section. SPM, SPMe and DFN do not merely differ in accuracy — they express different physics. Data generated from the wrong model has labels with no corresponding physical meaning.
- SPM: assumes infinite electrolyte conductivity and uniform concentration. Fine for macroscopic SOH or simple RUL when electrolyte dynamics are not limiting. At high C-rate it is systematically optimistic — a bias, not random error.
- SPMe: reintroduces an analytical approximation of electrolyte concentration gradients. Captures the dominant first-order effect at a fraction of DFN’s cost, which makes it the practical choice for broad parameter sweeps.
- DFN: resolves spatial distribution across anode, separator and cathode. Irreplaceable when the target is high-rate polarisation, localised plating, or EIS arcs where solid-electrolyte coupling dominates.
A workable combination: sweep with SPMe, re-run the selected subset with DFN. This only holds together if the model used is recorded in metadata — otherwise the two batches merge and can never be separated again.
6. Schema design for AI samples
Design the dataset for falsifiability. An auditable sample carries:
- Boundary and initial conditions: initial $c_s$ distribution (SOC), thermal field, rate limits, frequency excitation grid.
- Observable features: terminal voltage statistics, $Z_{re}$ and $Z_{im}$, high-frequency intercept, Warburg tail coefficient.
- Embedded physical labels: moles of lithium consumed by SEI (LLI), active material volume fraction lost to particle fracture (LAM, per electrode), local ECM parameters.
- Solver fingerprint: solver type, tolerance, version, failure flags, step statistics.
7. Sample audit table
Use this to judge whether a sample has research value, rather than merely producing a curve.
| Audit dimension | Must record | Research value | Unacceptable signal |
|---|---|---|---|
| Model structure | SPM/SPMe/DFN, enabled degradation submodels, parameter set version | Determines which mechanisms the sample can express | Only a voltage curve, with no record of model or parameter set |
| Experimental protocol | C-rate, temperature, cutoffs, rest time, EIS frequency grid | Determines identifiability between observables and internal state | Train and test share adjacent cycles of one trajectory |
| Solver state | Solver type and version, tolerance, failure flags, step statistics | Separates real physical change from numerical error | Non-convergent, deeply degraded samples written into training |
| Trajectory completeness | Cycles completed vs cycles requested | Identifies truncated, early-terminated trajectories | Trajectories of different lengths treated as equivalent samples |
| Label extraction | Formulas and units for SOH, LLI, LAM, ECM parameters | Ensures the model learns interpretable physical quantities | A bare health_score with no unit, formula or provenance |
8. Three methodological traps
- Treating simulation as ground truth. DFN output is a projection of one mathematical theory. It has to be calibrated against real cycler data by nonlinear least squares.
- Chasing sample count over physical diversity. Naive Monte Carlo sampling produces 10,000 nearly identical curves and collapses the data manifold. Sample with Latin hypercube across orthogonal physical dimensions.
- Ignoring solver tolerance. Extracting impedance features from deeply degraded states with inappropriate tolerances injects numerical noise indistinguishable from electrochemistry.
And one more fundamental risk: the danger of synthetic data is not inaccuracy but self-consistency. Real cells have manufacturing spread, thermal non-uniformity and sensor noise; a PDE solution has none. Metrics that look excellent in simulation can collapse outright on hardware.
References
- Python Battery Mathematical Modelling (PyBaMM), Journal of Open Research Software
- PyBaMM EIS Simulation documentation
- PyBaMM coupled degradation notebook
- pybamm-eis (archived; merged into core)
The “solver dies at cycle 300 in the middle of the night” scenario from the opening has its own troubleshooting guide: PyBaMM solver convergence failures – numerics or physics, including every IDAKLU option default and one trap that fails silently.