From d4ac1c3636eba8973f80217150b4179b6e874e0d Mon Sep 17 00:00:00 2001 From: David Straub Date: Wed, 30 Sep 2026 16:54:28 +0200 Subject: [PATCH] Solve DAE models with PathSim's solvers in CellElectrical and CellElectrothermal --- README.md | 39 ++--- pyproject.toml | 2 +- src/pathsim_batt/cells/pybamm_cell.py | 218 +++++++++++++++----------- tests/cells/test_lead_acid.py | 67 +++----- tests/cells/test_pybamm_cell.py | 85 ++++++++-- 5 files changed, 229 insertions(+), 182 deletions(-) diff --git a/README.md b/README.md index 0ae2f3d..6fb7a24 100644 --- a/README.md +++ b/README.md @@ -57,35 +57,25 @@ Two decisions determine the right block: **thermal ownership** and **integration | Block | Thermal | Strategy | Use when | |---|---|---|---| -| `CellElectrothermal` | PyBaMM (internal) | Monolithic ODE | Single cell, coupled electro-thermal, ODE model | -| `CellElectrical` + `LumpedThermal` | PathSim (external) | Monolithic ODE | Pack-level, custom cooling, ODE model | -| `CellCoSimElectrothermal` | PyBaMM (internal) | Co-simulation | DAE models (DFN, lead_acid.Full), mixed solvers | -| `CellCoSimElectrical` + `LumpedThermal` | PathSim (external) | Co-simulation | DAE models with external thermal network | +| `CellElectrothermal` | PyBaMM (internal) | Monolithic | Single cell, coupled electro-thermal | +| `CellElectrical` + `LumpedThermal` | PathSim (external) | Monolithic | Pack-level, custom cooling | +| `CellCoSimElectrothermal` | PyBaMM (internal) | Co-simulation | PyBaMM's own solvers, mixed solvers | +| `CellCoSimElectrical` + `LumpedThermal` | PathSim (external) | Co-simulation | PyBaMM's own solvers with external thermal network | `LumpedThermal` is a single-node thermal block (`mass`, `Cp`, `UA`, `T0`) that receives `Q_dot` from a `CellElectrical` block and feeds back cell temperature. -## PyBaMM model compatibility +## PyBaMM models -Thermal sub-model and heat-source options are injected automatically — pass the bare model class with no `options=`. +All blocks accept any PyBaMM battery model, e.g. `lithium_ion.SPM`, `SPMe`, `DFN`, `lead_acid.LOQS`, `lead_acid.Full` or `equivalent_circuit.Thevenin`. Thermal sub-model and heat-source options are injected automatically — pass the bare model class with no `options=`. -| PyBaMM model | Default parameter set | `CellElectrical` | `CellElectrothermal` | `CellCoSimElectrical` | `CellCoSimElectrothermal` | -|---|---|:---:|:---:|:---:|:---:| -| `lithium_ion.SPM` | `Chen2020` | ✅ | ✅ | ✅ | ✅ | -| `lithium_ion.SPMe` | `Chen2020` | ✅ | ✅ | ✅ | ✅ | -| `lithium_ion.DFN` | `Chen2020` | ❌ DAE | ❌ DAE | ✅ | ✅ | -| `lead_acid.LOQS` | `Sulzer2019` | ✅ ¹ | ✅ ¹ | ✅ ² | ✅ ² | -| `lead_acid.Full` | `Sulzer2019` | ❌ DAE | ❌ DAE | ✅ | ✅ | -| `equivalent_circuit.Thevenin` | `ECM_Example` | ✅ | ✅ | ✅ ³ | ✅ ³ | +Known limitations of the co-simulation blocks: -¹ Not on PyBaMM 26.7 — there `LOQS` is a DAE, use a `CellCoSim*` block instead. It is an ODE again from 26.8 on. - -² PyBaMM < 26.7 only — pass `pybamm_solver=pybamm.CasadiSolver(mode="safe")`; the default `IDAKLUSolver` errors on `LOQS`. Fixed in 26.7. - -³ `initial_soc=1.0` fails because PyBaMM requires event values to be strictly positive at `t=0`; the "Maximum SoC" event is zero exactly at full charge. Any value below 1.0 (e.g. `initial_soc=0.99`) works. +- `lead_acid.LOQS` on PyBaMM < 26.7: pass `pybamm_solver=pybamm.CasadiSolver(mode="safe")`; the default `IDAKLUSolver` errors on `LOQS`. +- `equivalent_circuit.Thevenin`: `initial_soc=1.0` fails because PyBaMM requires event values to be strictly positive at `t=0`; use e.g. `initial_soc=0.99`. ```python import pybamm -from pathsim_batt import CellElectrothermal, CellCoSimElectrical +from pathsim_batt import CellElectrical, CellElectrothermal # Custom chemistry / parameter set cell = CellElectrothermal( @@ -93,11 +83,10 @@ cell = CellElectrothermal( parameter_values=pybamm.ParameterValues("Mohtat2020"), ) -# Lead-acid via co-simulation (DAE model) -cell = CellCoSimElectrical( - model=pybamm.lead_acid.Full(), - parameter_values=pybamm.ParameterValues("Sulzer2019"), - dt=1.0, +# High-fidelity model +cell = CellElectrical( + model=pybamm.lithium_ion.DFN(), + parameter_values=pybamm.ParameterValues("Chen2020"), ) # Equivalent circuit model diff --git a/pyproject.toml b/pyproject.toml index 31eee7e..03ed180 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,7 @@ classifiers = [ "Topic :: Scientific/Engineering", ] dependencies = [ - "pathsim>=0.22", + "pathsim>=0.26", # PyBaMM ships native code that can't be installed in Pyodide. The # marker resolves to a normal eager install on every real platform and # is skipped on Emscripten so the pure-Python parts of the toolbox diff --git a/src/pathsim_batt/cells/pybamm_cell.py b/src/pathsim_batt/cells/pybamm_cell.py index 5f2c296..c5eee09 100644 --- a/src/pathsim_batt/cells/pybamm_cell.py +++ b/src/pathsim_batt/cells/pybamm_cell.py @@ -13,8 +13,9 @@ import numpy as np import numpy.typing as npt import pybamm -from pathsim.blocks import DynamicalSystem, Wrapper +from pathsim.blocks import Block, SemiExplicitDAE, Wrapper from pathsim.exceptions import StopSimulation +from pathsim.utils.register import Register # HELPERS ============================================================================= @@ -235,22 +236,21 @@ def _build_simulation( # BLOCKS =============================================================================== -class _CellBase(DynamicalSystem): +class _CellBase(SemiExplicitDAE): """Shared base for PyBaMM cell blocks. - Discretises the PyBaMM model at construction time and exposes its ODE - right-hand side to PathSim's numerical integrator via the ``DynamicalSystem`` - interface. The differential state vector of the discretised model becomes - the PathSim state; PathSim's chosen solver advances it in time. + Discretises the PyBaMM model at construction time and exposes it to + PathSim's numerical integrator via the ``SemiExplicitDAE`` interface. The + differential state vector of the discretised model becomes the PathSim + state; PathSim's chosen solver advances it in time. Algebraic variables + of models that result in a DAE system (e.g. DFN, lead_acid.Full) are + eliminated at every evaluation by solving the algebraic equations with + analytic CasADi Jacobians. Pure ODE models (e.g. SPMe, SPM) have no + algebraic variables. - Only PyBaMM models that produce a pure ODE after discretisation are - supported (i.e. models with no algebraic variables, such as SPMe and SPM). - Models that result in a DAE system (e.g. DFN) are not supported and will - raise ``NotImplementedError`` at construction time. - - Because the SPMe/SPM family of models is stiff, users should prefer an - implicit solver (e.g. ``ESDIRK43``, ``BDF``) when constructing the - PathSim ``Simulation``. + Because the discretised models are stiff, users should prefer an implicit + solver (e.g. ``ESDIRK43``, ``BDF``) when constructing the PathSim + ``Simulation``. Subclasses set ``_thermal_option`` and ``_pybamm_output_vars`` to select the thermal sub-model and define which PyBaMM variables map to the block's @@ -271,6 +271,7 @@ def __init__( parameter_values: pybamm.ParameterValues | None = None, initial_soc: float = 1.0, pybamm_solver: pybamm.BaseSolver | None = None, + tolerance: float = 1e-6, ) -> None: self._initial_soc = float(initial_soc) @@ -313,22 +314,6 @@ def __init__( available = sim.built_model.variables - # Early DAE check: probe with just the voltage variable so that the - # NotImplementedError is raised before variable-resolution, giving a - # cleaner error message for models that are DAE *and* also lack other - # required output variables (e.g. sodium_ion.BasicDFN). - _vol_var = _pick_var(available, _VOLTAGE_VAR_CANDIDATES, "terminal voltage") - _probe = sim.built_model.export_casadi_objects( - [_vol_var], input_parameter_order=list(_DEFAULT_INPUTS.keys()) - ) - if _probe["z"].numel() > 0: - raise NotImplementedError( - f"{type(self).__name__}: the supplied PyBaMM model has " - f"{_probe['z'].numel()} algebraic variable(s) after discretisation " - "(DAE system). Only pure ODE models are supported by this block. " - "Use a CellCoSim* block for DAE models." - ) - resolved_output_vars, soc_cap_var, soc_direct_var = _resolve_output_vars( self._pybamm_output_vars, available ) @@ -343,82 +328,109 @@ def __init__( input_parameter_order=list(_DEFAULT_INPUTS.keys()), ) - t_sym = casadi_objs["t"] - x_sym = casadi_objs["x"] - p_sym = casadi_objs["inputs"] - - rhs_fn = casadi.Function("rhs", [t_sym, x_sym, p_sym], [casadi_objs["rhs"]]) - jac_fn = casadi.Function( - "jac_rhs", [t_sym, x_sym, p_sym], [casadi_objs["jac_rhs"]] + args = [ + casadi_objs["t"], + casadi_objs["x"], + casadi_objs["z"], + casadi_objs["inputs"], + ] + n_x = casadi_objs["x"].numel() + jac_rhs = casadi_objs["jac_rhs"] + jac_alg = casadi_objs["jac_algebraic"] + + rhs_fn = casadi.Function("rhs", args, [casadi_objs["rhs"]]) + alg_fn = casadi.Function("alg", args, [casadi_objs["algebraic"]]) + jac_fns = [ + casadi.Function(name, args, [expr]) + for name, expr in [ + ("jac_rhs_x", jac_rhs[:, :n_x]), + ("jac_rhs_z", jac_rhs[:, n_x:]), + ("jac_alg_x", jac_alg[:, :n_x]), + ("jac_alg_z", jac_alg[:, n_x:]), + ] + ] + out_fn = casadi.Function( + "out", args, [casadi_objs["variables"][n] for n in all_out_vars] ) - out_var_fns = {} - for idx, var_name in enumerate(all_out_vars): - var_expr = casadi_objs["variables"][var_name] - out_var_fns[var_name] = casadi.Function( - f"outvar_{idx}", [t_sym, x_sym, p_sym], [var_expr] - ) - - self._casadi_rhs = rhs_fn - self._jac_rhs_eval = jac_fn - self._out_var_fcns = out_var_fns + self._resolved_output_vars = resolved_output_vars + self._soc_cap_var = soc_cap_var + self._out_fn = out_fn self._q_nominal = float(self._parameter_values["Nominal cell capacity [A.h]"]) - - q_nominal = self._q_nominal - initial_soc_val = float(initial_soc) - soc_direct_scale = _detect_soc_direct_scale(sim, soc_direct_var) + self._soc_direct_scale = _detect_soc_direct_scale(sim, soc_direct_var) def _pack(u): return casadi.DM([float(u[0]), float(u[1])]) - def func_dyn(x, u, t): - xv = casadi.DM(x.reshape(-1, 1)) - p = _pack(u) - return np.array(rhs_fn(t, xv, p)).flatten() - - def jac_dyn(x, u, t): - xv = casadi.DM(x.reshape(-1, 1)) - p = _pack(u) - return np.array(jac_fn(t, xv, p)) - - v_lower = self._v_lower - v_upper = self._v_upper - v_idx = self._v_idx - - def func_alg(x, u, t): - xv = casadi.DM(x.reshape(-1, 1)) - p = _pack(u) - outputs = [float(out_var_fns[n](t, xv, p)) for n in resolved_output_vars] - if soc_cap_var is not None: - q_dis = float(out_var_fns[soc_cap_var](t, xv, p)) - soc = max(0.0, min(1.0, initial_soc_val - q_dis / q_nominal)) - else: - raw = float(out_var_fns[soc_direct_var](t, xv, p)) - soc = max(0.0, min(1.0, raw * soc_direct_scale)) - outputs.append(soc) - V = outputs[v_idx] - if V <= v_lower: - raise StopSimulation(f"undervoltage: V={V:.4f} V <= {v_lower} V") - if V >= v_upper: - raise StopSimulation(f"overvoltage: V={V:.4f} V >= {v_upper} V") - return np.array(outputs) - - x0_fn = casadi.Function("x0", [p_sym], [casadi_objs["x0"]]) - - y0 = np.array(x0_fn(casadi.DM(list(_DEFAULT_INPUTS.values())))).flatten() + def func_dyn(x, z, u, t): + return np.array(rhs_fn(t, x, z, _pack(u))).ravel() + + def func_alg(x, z, u, t): + return np.array(alg_fn(t, x, z, _pack(u))).ravel() + + # CasADi's sparse export is much faster than a dense DM conversion + # for the large Jacobians of spatially resolved models. + def _jac(fn): + return lambda x, z, u, t: fn(t, x, z, _pack(u)).sparse().toarray() + + p0 = casadi.DM(list(_DEFAULT_INPUTS.values())) + x0 = casadi.Function("x0", [args[3]], [casadi_objs["x0"]])(p0) + z0 = casadi.Function("z0", [args[3]], [casadi_objs["z0"]])(p0) super().__init__( func_dyn=func_dyn, func_alg=func_alg, - initial_value=y0, - jac_dyn=jac_dyn, + initial_value=np.array(x0).ravel(), + z0=np.array(z0).ravel(), + jac_dyn_x=_jac(jac_fns[0]), + jac_dyn_z=_jac(jac_fns[1]), + jac_alg_x=_jac(jac_fns[2]), + jac_z=_jac(jac_fns[3]), + tolerance=tolerance, ) + # SemiExplicitDAE sizes its outputs to the stacked state [x, z]; this + # block exposes the cell outputs instead. + self.outputs = Register( + size=len(self.output_port_labels), + mapping=self.output_port_labels.copy(), + ) + + def _cell_outputs(self, x, z, u, t) -> npt.NDArray[np.float64]: + """Evaluate the output variables and the SOC at the given state.""" + values = [float(v) for v in self._out_fn(t, x, z, casadi.DM(u[:2]))] + outputs = values[:-1] + if self._soc_cap_var is not None: + soc = self._initial_soc - values[-1] / self._q_nominal + else: + soc = values[-1] * self._soc_direct_scale + outputs.append(max(0.0, min(1.0, soc))) + return np.array(outputs) + + def update(self, t: float) -> None: + """Eliminate the algebraic states and evaluate the cell outputs. + + Raises ``StopSimulation`` when the terminal voltage leaves the + cut-off window of the parameter set. + """ + x, u = self.engine.state, self.inputs.to_array() + self._z = self._solve_z(x, u, t) + outputs = self._cell_outputs(x, self._z, u, t) + self.outputs.update_from_array(outputs) + V = outputs[self._v_idx] + if V <= self._v_lower: + raise StopSimulation(f"undervoltage: V={V:.4f} V <= {self._v_lower} V") + if V >= self._v_upper: + raise StopSimulation(f"overvoltage: V={V:.4f} V >= {self._v_upper} V") + def __len__(self) -> int: return len(self._pybamm_output_vars) + 1 def reset(self) -> None: - super().reset() + # Bypass SemiExplicitDAE.reset, which writes the stacked state + # [x, z] to the outputs. + Block.reset(self) + self._z = self.z0.copy() class _CoSimCellBase(Wrapper): @@ -607,12 +619,16 @@ class CellElectrical(_CellBase): """Cell block — electrical outputs only, external thermal coupling. PathSim integrates the electrochemical state via the discretised PyBaMM - ODE. Temperature dynamics live outside this block: wire ``Q_dot`` to a + model. Temperature dynamics live outside this block: wire ``Q_dot`` to a ``LumpedThermal`` (or similar) block and feed its temperature output back to ``T_cell``. + Models that result in a DAE system after discretisation (e.g. DFN, + lead_acid.Full) are supported; their algebraic variables are solved by + the block at every evaluation. + .. note:: - The SPMe/SPM ODE is stiff. Use an implicit solver (e.g. + The discretised models are stiff. Use an implicit solver (e.g. ``ESDIRK43``, ``BDF``) when constructing the PathSim ``Simulation`` to avoid prohibitively small step sizes. @@ -628,6 +644,9 @@ class CellElectrical(_CellBase): pybamm_solver : pybamm.BaseSolver or None PyBaMM solver used only during model build / discretisation. Defaults to ``IDAKLUSolver()``. + tolerance : float + Convergence tolerance on the residual norm of the algebraic equations + of DAE models. Default 1e-6. Inputs ------ @@ -656,13 +675,17 @@ class CellElectrothermal(_CellBase): """Cell block — coupled electrical and thermal model. PathSim integrates the full electrochemical + thermal state (via the - discretised PyBaMM ODE). The cell temperature is part of the PyBaMM + discretised PyBaMM model). The cell temperature is part of the PyBaMM state vector and is read back as output port ``T``. Supply a time-varying ambient / coolant temperature via ``T_amb`` to couple to a pack-level thermal model. + Models that result in a DAE system after discretisation (e.g. DFN, + lead_acid.Full) are supported; their algebraic variables are solved by + the block at every evaluation. + .. note:: - The SPMe/SPM ODE is stiff. Use an implicit solver (e.g. + The discretised models are stiff. Use an implicit solver (e.g. ``ESDIRK43``, ``BDF``) when constructing the PathSim ``Simulation`` to avoid prohibitively small step sizes. @@ -677,6 +700,9 @@ class CellElectrothermal(_CellBase): pybamm_solver : pybamm.BaseSolver or None PyBaMM solver used only during model build / discretisation. Defaults to ``IDAKLUSolver()``. + tolerance : float + Convergence tolerance on the residual norm of the algebraic equations + of DAE models. Default 1e-6. Inputs ------ @@ -709,7 +735,8 @@ class CellCoSimElectrical(_CoSimCellBase): ``pybamm.Simulation.step()``. PathSim receives zero-order-held outputs between macro-steps. - This mode supports PyBaMM models that result in DAE systems (e.g. DFN). + PyBaMM's own solvers handle any PyBaMM model, including those that + result in DAE systems (e.g. DFN). Parameters ---------- @@ -745,7 +772,8 @@ class CellCoSimElectrothermal(_CoSimCellBase): ``pybamm.Simulation.step()``. PathSim receives zero-order-held outputs between macro-steps. - This mode supports PyBaMM models that result in DAE systems (e.g. DFN). + PyBaMM's own solvers handle any PyBaMM model, including those that + result in DAE systems (e.g. DFN). Parameters ---------- diff --git a/tests/cells/test_lead_acid.py b/tests/cells/test_lead_acid.py index 633143b..6ce8a39 100644 --- a/tests/cells/test_lead_acid.py +++ b/tests/cells/test_lead_acid.py @@ -2,9 +2,8 @@ Block / model matrix covered ----------------------------- -lead_acid.LOQS — ODE (all 4 blocks), except on PyBaMM 26.7 where it is a DAE - → CoSim blocks only there -lead_acid.Full — DAE → CoSim blocks only +lead_acid.LOQS — ODE (DAE on PyBaMM 26.7), all 4 blocks +lead_acid.Full — DAE, all 4 blocks """ import unittest @@ -18,7 +17,6 @@ CellCoSimElectrical, CellCoSimElectrothermal, CellElectrical, - CellElectrothermal, ) from ._helpers import ( @@ -30,24 +28,16 @@ run_electrothermal, ) -# PyBaMM 26.7 registers "voltage as a state" centrally on every -# BaseBatteryModel (default "true"), and lead-acid models don't support -# disabling it, so LOQS is a DAE there and can't run in the monolithic -# (ODE-only) blocks. 26.8 defaults it back to "false", making LOQS an ODE -# again. Detect this from the model itself rather than the version number. -_LOQS_IS_ODE = not pybamm.lead_acid.LOQS().algebraic - # --------------------------------------------------------------------------- -# lead_acid.LOQS (ODE — all 4 blocks) +# lead_acid.LOQS # --------------------------------------------------------------------------- class TestLeadAcidLOQS(unittest.TestCase): """lead_acid.LOQS with Sulzer2019 parameters. - ODE model (all 4 blocks), except on PyBaMM 26.7 where it is a DAE (CoSim - blocks only). Sulzer2019 cutoffs: lower 1.75 V, upper 2.42 V, nominal capacity - 17 A·h. + ODE model, except on PyBaMM 26.7 where it is a DAE. Sulzer2019 cutoffs: + lower 1.75 V, upper 2.42 V, nominal capacity 17 A·h. """ def setUp(self): @@ -58,28 +48,14 @@ def setUp(self): def _model(self): return pybamm.lead_acid.LOQS() - @unittest.skipUnless(_LOQS_IS_ODE, "LOQS is a DAE on this PyBaMM version") def test_electrical_smoke(self): cell = run_electrical(self._model(), self.pv, current=17.0) assert_electrical_outputs(self, cell, self.v_lo, self.v_hi) - @unittest.skipUnless(_LOQS_IS_ODE, "LOQS is a DAE on this PyBaMM version") def test_electrothermal_smoke(self): cell = run_electrothermal(self._model(), self.pv, current=17.0) assert_electrothermal_outputs(self, cell, self.v_lo, self.v_hi) - @unittest.skipIf(_LOQS_IS_ODE, "LOQS is a pure ODE on this PyBaMM version") - def test_monolithic_electrical_raises(self): - """Where LOQS is a DAE (PyBaMM 26.7), CellElectrical must raise.""" - with self.assertRaises(NotImplementedError): - CellElectrical(model=self._model(), parameter_values=self.pv) - - @unittest.skipIf(_LOQS_IS_ODE, "LOQS is a pure ODE on this PyBaMM version") - def test_monolithic_electrothermal_raises(self): - """Where LOQS is a DAE (PyBaMM 26.7), CellElectrothermal must raise.""" - with self.assertRaises(NotImplementedError): - CellElectrothermal(model=self._model(), parameter_values=self.pv) - def test_cosim_electrical_smoke(self): # On PyBaMM < 26.7, LOQS disables its Jacobian and IDAKLUSolver (the # co-sim default) errors without one; CasadiSolver works on every @@ -127,19 +103,16 @@ def test_cosim_electrothermal_smoke(self): sim.run(2) assert_electrothermal_outputs(self, cell, self.v_lo, self.v_hi) - @unittest.skipUnless(_LOQS_IS_ODE, "LOQS is a DAE on this PyBaMM version") def test_electrical_soc_decreases(self): """SOC must decrease under discharge current.""" cell = run_electrical(self._model(), self.pv, current=17.0, duration=60) self.assertLess(float(cell.outputs[2]), 1.0) - @unittest.skipUnless(_LOQS_IS_ODE, "LOQS is a DAE on this PyBaMM version") def test_cutoff_values_match_parameter_set(self): cell = CellElectrical(model=self._model(), parameter_values=self.pv) self.assertAlmostEqual(cell._v_lower, self.v_lo) self.assertAlmostEqual(cell._v_upper, self.v_hi) - @unittest.skipUnless(_LOQS_IS_ODE, "LOQS is a DAE on this PyBaMM version") def test_q_dot_nonzero_during_discharge(self): """Q_dot must be strictly positive during discharge (isothermal LOQS). @@ -153,7 +126,6 @@ def test_q_dot_nonzero_during_discharge(self): "Q_dot is zero — thermal model may not compute heat sources", ) - @unittest.skipUnless(_LOQS_IS_ODE, "LOQS is a DAE on this PyBaMM version") def test_tamb_affects_temperature(self): """A warmer ambient temperature must yield a higher output cell temperature.""" solver = pybamm.CasadiSolver(mode="safe") @@ -184,7 +156,6 @@ def test_tamb_affects_temperature(self): ), ) - @unittest.skipUnless(_LOQS_IS_ODE, "LOQS is a DAE on this PyBaMM version") def test_soc_scale_factor(self): """SOC must be well below 1.0 after sustained discharge. @@ -200,12 +171,12 @@ def test_soc_scale_factor(self): # --------------------------------------------------------------------------- -# lead_acid.Full (DAE — co-simulation only) +# lead_acid.Full (DAE) # --------------------------------------------------------------------------- class TestLeadAcidFull(unittest.TestCase): - """lead_acid.Full with Sulzer2019 parameters (DAE model — co-sim only).""" + """lead_acid.Full with Sulzer2019 parameters (DAE model).""" def setUp(self): self.pv = pybamm.ParameterValues("Sulzer2019") @@ -215,15 +186,23 @@ def setUp(self): def _model(self): return pybamm.lead_acid.Full() - def test_monolithic_electrical_raises(self): - """Full is a DAE — CellElectrical must raise NotImplementedError.""" - with self.assertRaises(NotImplementedError): - CellElectrical(model=self._model(), parameter_values=self.pv) + def test_electrical_smoke(self): + cell = run_electrical(self._model(), self.pv, current=17.0) + assert_electrical_outputs(self, cell, self.v_lo, self.v_hi) + self.assertGreater(len(cell.z0), 0) + + def test_electrothermal_smoke(self): + cell = run_electrothermal(self._model(), self.pv, current=17.0) + assert_electrothermal_outputs(self, cell, self.v_lo, self.v_hi) - def test_monolithic_electrothermal_raises(self): - """Full is a DAE — CellElectrothermal must raise NotImplementedError.""" - with self.assertRaises(NotImplementedError): - CellElectrothermal(model=self._model(), parameter_values=self.pv) + def test_electrical_matches_pybamm(self): + """Voltage after 10 min at 1 A must match PyBaMM's own solver.""" + cell = run_electrical(self._model(), self.pv, current=1.0, duration=600) + pv = self.pv.copy() + pv["Current function [A]"] = 1.0 + sol = pybamm.Simulation(self._model(), parameter_values=pv).solve([0, 600]) + V_ref = float(sol["Voltage [V]"].entries[-1]) + self.assertAlmostEqual(float(cell.outputs[0]), V_ref, delta=1e-3) def test_cosim_electrical_smoke(self): cell = run_cosim_electrical(self._model(), self.pv, current=17.0) diff --git a/tests/cells/test_pybamm_cell.py b/tests/cells/test_pybamm_cell.py index 0db97cf..0715e15 100644 --- a/tests/cells/test_pybamm_cell.py +++ b/tests/cells/test_pybamm_cell.py @@ -13,6 +13,8 @@ CellElectrothermal, ) +from ._helpers import assert_electrothermal_outputs, run_electrical, run_electrothermal + class TestPorts(unittest.TestCase): def test_electrical_input_labels(self): @@ -71,11 +73,13 @@ def test_initial_value_is_numpy_array(self): self.assertIsInstance(cell.initial_value, np.ndarray) self.assertGreater(len(cell.initial_value), 1) - def test_has_casadi_rhs(self): - """CasADi RHS is compiled and callable at construction time.""" + def test_rhs_matches_state_size(self): + """The reduced right-hand side returns one derivative per state.""" for cls in (CellElectrical, CellElectrothermal): cell = cls() - self.assertIsNotNone(cell._casadi_rhs) + u = np.array([0.0, 298.15]) + dx = cell.op_dyn(cell.initial_value, u, 0.0) + self.assertEqual(dx.shape, cell.initial_value.shape) def test_state_size_equals_differential_states_only(self): """State must contain only differential (x) variables, not algebraic (z).""" @@ -108,27 +112,43 @@ def test_state_size_equals_differential_states_only(self): expected_x_size = objs["x"].numel() self.assertEqual(len(cell.initial_value), expected_x_size) - def test_jac_dyn_is_square(self): - """jac_dyn must return a square (n×n) matrix where n is the state size.""" + def test_jac_is_square(self): + """The reduced Jacobian must be a square (n×n) matrix, n the state size.""" for cls in (CellElectrical, CellElectrothermal): cell = cls() n = len(cell.initial_value) - x = cell.initial_value u = np.array([0.0, 298.15]) - J = cell.jac_dyn(x, u, 0.0) + J = cell.op_dyn.jac_x(cell.initial_value, u, 0.0) self.assertEqual(J.shape, (n, n)) - def test_dfn_model_raises(self): - """DFN models (DAE after discretisation) must raise NotImplementedError.""" + def test_spme_has_no_algebraic_states(self): + self.assertEqual(len(CellElectrical().z0), 0) + + def test_dfn_supported(self): + """DFN (DAE after discretisation) has algebraic states and is supported.""" dfn = pybamm.lithium_ion.DFN(options={"thermal": "isothermal"}) - with self.assertRaises(NotImplementedError): - CellElectrical(model=dfn) + cell = CellElectrical(model=dfn) + self.assertGreater(len(cell.z0), 0) + self.assertEqual(len(cell), 3) - def test_dfn_lumped_raises(self): - """DFN with lumped thermal also has algebraic variables and must raise.""" + def test_dfn_lumped_supported(self): + """DFN with lumped thermal is supported by the electrothermal block.""" dfn = pybamm.lithium_ion.DFN(options={"thermal": "lumped"}) - with self.assertRaises(NotImplementedError): - CellElectrothermal(model=dfn) + cell = CellElectrothermal(model=dfn) + self.assertGreater(len(cell.z0), 0) + self.assertEqual(len(cell), 4) + + def test_outputs_sized_to_ports(self): + """Outputs are the cell ports, not the stacked state [x, z].""" + dfn = pybamm.lithium_ion.DFN(options={"thermal": "isothermal"}) + cell = CellElectrical(model=dfn) + self.assertEqual(len(cell.outputs), 3) + cell.reset() + self.assertEqual(len(cell.outputs), 3) + + def test_tolerance_passed_through(self): + self.assertEqual(CellElectrical().tolerance, 1e-6) + self.assertEqual(CellElectrical(tolerance=1e-8).tolerance, 1e-8) def test_dfn_cosim_supported(self): """DFN is supported by co-simulation blocks.""" @@ -144,7 +164,7 @@ def test_dfn_cosim_electrothermal_supported(self): class TestElectrical(unittest.TestCase): - """Integration tests for CellElectrical — PathSim integrates the PyBaMM ODE.""" + """Integration tests for CellElectrical — PathSim integrates the PyBaMM model.""" def _make_simulation(self, cell, current, T_cell): """Create a Simulation with the cell and constant inputs.""" @@ -396,6 +416,37 @@ def _run_and_get_T_cell(T_amb): ) +class TestDFN(unittest.TestCase): + """DFN (DAE) integrated by PathSim, compared against PyBaMM's own solver.""" + + def _reference_voltage(self, model, current, t_end): + pv = pybamm.ParameterValues("Chen2020") + pv["Current function [A]"] = current + sim = pybamm.Simulation(model, parameter_values=pv) + sol = sim.solve([0, t_end], initial_soc=1.0) + return float(sol["Voltage [V]"].entries[-1]) + + def test_electrical_matches_pybamm(self): + dfn = pybamm.lithium_ion.DFN() + cell = run_electrical(dfn, pybamm.ParameterValues("Chen2020"), 5.0, 298.15, 60) + V_ref = self._reference_voltage(dfn, 5.0, 60) + self.assertAlmostEqual(float(cell.outputs[0]), V_ref, delta=1e-3) + self.assertLess(float(cell.outputs[2]), 1.0) + self.assertGreater(float(cell.outputs[1]), 0.0) + + def test_electrothermal_outputs_physical(self): + dfn = pybamm.lithium_ion.DFN() + pv = pybamm.ParameterValues("Chen2020") + cell = run_electrothermal(dfn, pv, 5.0, 298.15, 60) + assert_electrothermal_outputs( + self, + cell, + float(pv["Lower voltage cut-off [V]"]), + float(pv["Upper voltage cut-off [V]"]), + ) + self.assertLess(float(cell.outputs[3]), 1.0) + + class TestCoSimulationElectrical(unittest.TestCase): """Integration tests for CellCoSimElectrical — PyBaMM performs the stepping.""" @@ -656,7 +707,7 @@ def _run(self, cell, current, T_input_port, T_value, cosim_dt=None): def test_non_cosim_stops_before_negative_voltage(self): """CellElectrical must stop automatically before V goes negative. - ``StopSimulation`` is raised from ``func_alg`` the moment voltage + ``StopSimulation`` is raised from ``update`` the moment voltage reaches the lower cut-off, so PathSim halts without any user wiring. """ cell = CellElectrical(initial_soc=0.02)