diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 3a75104..fdf1f08 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -27,8 +27,18 @@ jobs: test: runs-on: ubuntu-latest strategy: + fail-fast: false matrix: - python-version: ["3.10", "3.11", "3.12", "3.13"] + include: + # All Python versions against the latest PyBaMM + - {python-version: "3.10", pybamm-version: "26.4"} + - {python-version: "3.11", pybamm-version: "26.4"} + - {python-version: "3.12", pybamm-version: "26.4"} + - {python-version: "3.13", pybamm-version: "26.4"} + - {python-version: "3.14", pybamm-version: "26.4"} + # Older PyBaMM versions against the latest Python they support + - {python-version: "3.13", pybamm-version: "25.12"} + - {python-version: "3.14", pybamm-version: "26.3"} steps: - uses: actions/checkout@v4 @@ -38,7 +48,7 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install - run: pip install -e ".[test]" mypy + run: pip install -e ".[test]" mypy "pybamm==${{ matrix.pybamm-version }}" - name: mypy run: mypy src/pathsim_batt diff --git a/README.md b/README.md index a229839..9f83033 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ --- -PathSim-Batt extends the [PathSim](https://github.com/pathsim/pathsim) simulation framework with battery cell blocks using [PyBaMM](https://pybamm.org) as the electrochemical backend. All blocks follow the standard PathSim block interface and can be connected into simulation diagrams. +PathSim-Batt extends the [PathSim](https://github.com/pathsim/pathsim) simulation framework with battery cell blocks backed by [PyBaMM](https://pybamm.org). All blocks follow the standard PathSim interface and can be wired into any simulation diagram. ## Install @@ -28,54 +28,83 @@ PathSim-Batt extends the [PathSim](https://github.com/pathsim/pathsim) simulatio pip install pathsim-batt ``` -## Blocks +## Quick start -| Block | Description | Key Parameters | -|-------|-------------|----------------| -| `CellElectrothermal` | Coupled electrical + thermal cell (PathSim integrates PyBaMM ODE incl. temperature) | `model`, `parameter_values`, `initial_soc` | -| `CellElectrical` | Electrical only, isothermal; wire to `LumpedThermal` for external thermal coupling | `model`, `parameter_values`, `initial_soc` | -| `CellCoSimElectrothermal` | Coupled electrical + thermal co-simulation cell (PyBaMM steps internally) | `model`, `parameter_values`, `initial_soc`, `dt` | -| `CellCoSimElectrical` | Electrical co-simulation cell for external thermal coupling | `model`, `parameter_values`, `initial_soc`, `dt` | -| `LumpedThermal` | Single-node thermal model for external thermal coupling | `mass`, `Cp`, `UA`, `T0` | +```python +import pybamm +from pathsim import Connection, Simulation +from pathsim.blocks import Constant +from pathsim.solvers import ESDIRK43 +from pathsim_batt import CellElectrothermal + +cell = CellElectrothermal(initial_soc=1.0) # defaults: SPMe + Chen2020 +I_src = Constant(5.0) # 5 A discharge +T_src = Constant(298.15) # 25 °C ambient + +sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[Connection(I_src, cell["I"]), Connection(T_src, cell["T_amb"])], + dt=1.0, + Solver=ESDIRK43, +) +sim.run(3600) +print(f"V = {cell.outputs[0]:.3f} V T = {cell.outputs[1]:.1f} K SOC = {cell.outputs[3]:.3f}") +``` + +## Choosing a block + +Two decisions determine the right block: **thermal ownership** and **integration strategy**. + +| 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 | -## PyBaMM integration +`LumpedThermal` is a single-node thermal block (`mass`, `Cp`, `UA`, `T0`) that receives `Q_dot` from a `CellElectrical` block and feeds back cell temperature. -The cell blocks wrap [PyBaMM](https://pybamm.org) models behind the PathSim block interface. +## PyBaMM model compatibility -- `CellElectrothermal` / `CellElectrical` use PathSim monolithic integration (`DynamicalSystem`) and exported CasADi ODE right-hand sides. -- `CellCoSimElectrothermal` / `CellCoSimElectrical` use periodic co-simulation (`Wrapper`) and call `pybamm.Simulation.step()` internally. +Thermal sub-model and heat-source options are injected automatically — pass the bare model class with no `options=`. -Only models that yield a **pure ODE** after discretisation are supported by the monolithic blocks (`CellElectrothermal`, `CellElectrical`) — currently SPMe and SPM. Models such as DFN that produce a DAE system (algebraic variables) will raise `NotImplementedError` there. +| 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` | ✅ | ✅ | ✅ ² | ✅ ² | -For DAE models (e.g. DFN), use the co-simulation blocks (`CellCoSimElectrothermal`, `CellCoSimElectrical`). +¹ Pass `pybamm_solver=pybamm.CasadiSolver(mode="safe")` — the default `IDAKLUSolver` requires a Jacobian that ODE models do not provide. -- **ODE-type PyBaMM models** (SPMe, SPM) can be injected via the `model` parameter -- **Any parameter set** can be used via `parameter_values` (defaults to `Chen2020`) -- **Immediate initialisation** — the PyBaMM model is discretised during block construction +² `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. ```python import pybamm +from pathsim_batt import CellElectrothermal, CellCoSimElectrical -model = pybamm.lithium_ion.SPMe(options={"thermal": "lumped"}) -params = pybamm.ParameterValues("Mohtat2020") -cell = CellElectrothermal(model=model, parameter_values=params) - -# DAE example (DFN): use co-simulation mode -dfn_cell = CellCoSimElectrothermal( - model=pybamm.lithium_ion.DFN(options={"thermal": "lumped"}), - parameter_values=params, - dt=0.1, +# Custom chemistry / parameter set +cell = CellElectrothermal( + model=pybamm.lithium_ion.SPMe(), + parameter_values=pybamm.ParameterValues("Mohtat2020"), ) -``` -## Thermal coupling modes +# Lead-acid via co-simulation (DAE model) +cell = CellCoSimElectrical( + model=pybamm.lead_acid.Full(), + parameter_values=pybamm.ParameterValues("Sulzer2019"), + dt=1.0, +) -| Mode | Block | Owns cell temperature | Use when | -|---|---|---|---| -| Internal | `CellElectrothermal` | PyBaMM | Single-cell simulations, quick setup | -| External | `CellElectrical` + `LumpedThermal` | PathSim | Multi-cell packs, custom cooling models | -| Co-sim internal | `CellCoSimElectrothermal` | PyBaMM | DAE models (e.g. DFN), mixed-solver workflows | -| Co-sim external | `CellCoSimElectrical` + `LumpedThermal` | PathSim | DAE models with external thermal network | +# Equivalent circuit model +cell = CellElectrical( + model=pybamm.equivalent_circuit.Thevenin(), + parameter_values=pybamm.ParameterValues("ECM_Example"), + initial_soc=0.9, +) +``` ## License diff --git a/pyproject.toml b/pyproject.toml index 4b44d6e..c0765f6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,7 @@ classifiers = [ "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", "Topic :: Scientific/Engineering", ] dependencies = [ diff --git a/src/pathsim_batt/cells/pybamm_cell.py b/src/pathsim_batt/cells/pybamm_cell.py index 3ec88c4..7f6d516 100644 --- a/src/pathsim_batt/cells/pybamm_cell.py +++ b/src/pathsim_batt/cells/pybamm_cell.py @@ -23,10 +23,29 @@ "Ambient temperature [K]": 298.15, } -# The PyBaMM variable name used for voltage cut-off detection. -# Both base classes locate it by name in ``_pybamm_output_vars`` at construction -# time, so subclasses may place it at any position in the list. -_TERMINAL_VOLTAGE_VAR = "Terminal voltage [V]" +_VOLTAGE_VAR_CANDIDATES: list[str] = [ + "Terminal voltage [V]", + "Voltage [V]", + "Battery voltage [V]", +] +_HEATING_VAR_CANDIDATES: list[str] = [ + "Total heating [W]", + "Total heat generation [W]", +] +_TEMP_VAR_CANDIDATES: list[str] = [ + "X-averaged cell temperature [K]", + "Cell temperature [K]", +] +_SOC_CAPACITY_VAR = "Discharge capacity [A.h]" +_SOC_DIRECT_CANDIDATES: list[str] = ["State of charge", "State of Charge", "SoC"] + +# Map from the canonical names used in ``_pybamm_output_vars`` to their +# per-model-family fallback lists. +_VAR_ALIAS_MAP: dict[str, list[str]] = { + "Terminal voltage [V]": _VOLTAGE_VAR_CANDIDATES, + "Total heating [W]": _HEATING_VAR_CANDIDATES, + "X-averaged cell temperature [K]": _TEMP_VAR_CANDIDATES, +} def _prepare_parameter_values( @@ -42,6 +61,156 @@ def _prepare_parameter_values( return parameter_values +def _pick_var(variables: dict, candidates: list[str], description: str) -> str: + """Return the first name in *candidates* present in *variables*. + + Raises ``ValueError`` listing all tried names if none match. + """ + for name in candidates: + if name in variables: + return name + raise ValueError( + f"No {description} variable found in PyBaMM model. Tried: {candidates}" + ) + + +def _resolve_output_vars( + requested: list[str], + available: dict, +) -> tuple[list[str], str | None, str | None]: + """Map *requested* variable names to names actually present in *available*. + + Returns ``(resolved, soc_cap_var, soc_direct_var)``: + + * ``resolved`` — same length as *requested*; each entry is the actual + PyBaMM variable name to use (canonical or alias). + * ``soc_cap_var`` — ``"Discharge capacity [A.h]"`` when available (standard + electrochemical models), otherwise ``None``. + * ``soc_direct_var`` — a direct SoC variable name used when *soc_cap_var* + is ``None`` (ECM, some basic models). + """ + resolved = [ + _pick_var(available, _VAR_ALIAS_MAP[name], name) + if (name not in available and name in _VAR_ALIAS_MAP) + else name + for name in requested + ] + + if _SOC_CAPACITY_VAR in available: + return resolved, _SOC_CAPACITY_VAR, None + return ( + resolved, + None, + _pick_var(available, _SOC_DIRECT_CANDIDATES, "state of charge"), + ) + + +def _detect_soc_direct_scale( + sim: pybamm.Simulation, + soc_direct_var: str | None, +) -> float: + """Return 1/100 if *soc_direct_var* is in percentage form, else 1.0. + + Some models (e.g. lead-acid) export ``"State of Charge"`` in the range + 0–100 rather than 0–1. Evaluating the variable at the initial state with + zero current lets us detect this once at construction time. + """ + if soc_direct_var is None: + return 1.0 + objs = sim.built_model.export_casadi_objects( + [soc_direct_var], + input_parameter_order=list(_DEFAULT_INPUTS.keys()), + ) + p0 = casadi.DM(list(_DEFAULT_INPUTS.values())) + x0 = casadi.Function("x0", [objs["inputs"]], [objs["x0"]])(p0) + z0 = casadi.Function("z0", [objs["inputs"]], [objs["z0"]])(p0) + soc_fn = casadi.Function( + "soc", + [objs["t"], objs["x"], objs["z"], objs["inputs"]], + [objs["variables"][soc_direct_var]], + ) + raw = float(soc_fn(0.0, x0, z0, p0)) + return 1.0 / 100.0 if raw > 1.0 else 1.0 + + +def _inject_thermal_options( + model: pybamm.BaseBatteryModel, + required_options: dict[str, str], +) -> pybamm.BaseBatteryModel: + """Return *model* with *required_options* merged into its options if needed. + + Handles both the thermal sub-model selection (``"thermal": "isothermal"`` + or ``"lumped"``) and ancillary flags such as + ``"calculate heat source for isothermal models": "true"``. This means + users can pass a plain ``pybamm.lead_acid.LOQS()`` to + ``CellElectrothermal`` without having to specify ``thermal='lumped'`` + themselves — the block injects it automatically. + + Models that already carry the required options are returned unchanged. + Models that have no ``"thermal"`` key in their options at all (e.g. ECM, + which manages temperature internally) are skipped silently — injection + does not apply to them. Models that expose ``"thermal"`` but whose + constructor rejects the specific option values emit a ``UserWarning`` + and are returned unchanged. + """ + if not required_options: + return model + # Models that have no "thermal" key in their options (e.g. ECM) manage + # temperature through their own internal mechanism and do not use PyBaMM's + # thermal sub-model system. Injection is not applicable; skip silently. + if "thermal" in required_options and "thermal" not in model.options: + return model + if all(model.options.get(k) == v for k, v in required_options.items()): + return model + try: + return type(model)(options={**dict(model.options), **required_options}) + except (pybamm.OptionError, TypeError): + import warnings + + warnings.warn( + f"{type(model).__name__} does not support options {required_options}; " + "thermal behaviour may be incorrect.", + UserWarning, + stacklevel=4, + ) + return model + + +def _build_simulation( + sim: pybamm.Simulation, + model: pybamm.BaseBatteryModel, + initial_soc: float, +) -> None: + """Build *sim*, handling model-family-specific initialisation quirks. + + * **ECM** (``Thevenin``): initial SoC is set via ``set_initial_state`` on + the parameter values before build; ``initial_soc`` cannot be passed to + ``sim.build()`` directly. + * **Lead-acid** models: the eSOH solver used by the standard build path is + lithium-ion–specific and fails for lead-acid parameter sets; build + without ``initial_soc``. + * **All other models**: standard ``sim.build(initial_soc=...)``. + """ + if isinstance(model, pybamm.equivalent_circuit.Thevenin): + from pybamm.models.full_battery_models.equivalent_circuit import ( + set_initial_state, + ) + + set_initial_state(initial_soc, sim.parameter_values) + sim.build(inputs=_DEFAULT_INPUTS) + elif isinstance(model, pybamm.lead_acid.BaseModel): + if initial_soc != 1.0: + raise ValueError( + "initial_soc is not supported for lead-acid models: PyBaMM's " + "lead-acid parameter sets do not include the electrode OCP data " + "required to map a target SoC to initial stoichiometries. " + "The initial state is always determined by the parameter values." + ) + sim.build(inputs=_DEFAULT_INPUTS) + else: + sim.build(initial_soc=initial_soc, inputs=_DEFAULT_INPUTS) + + # BLOCKS =============================================================================== @@ -65,8 +234,10 @@ class _CellBase(DynamicalSystem): Subclasses set ``_thermal_option`` and ``_pybamm_output_vars`` to select the thermal sub-model and define which PyBaMM variables map to the block's output ports (SOC is always appended last). ``_pybamm_output_vars`` must - contain ``_TERMINAL_VOLTAGE_VAR``; its position in the list is found - dynamically. + contain a terminal-voltage entry — either the canonical + ``"Terminal voltage [V]"`` or any alias listed in + ``_VOLTAGE_VAR_CANDIDATES``; the actual name exported by the built model is + resolved automatically at construction time. """ _thermal_option: str = "" @@ -82,18 +253,21 @@ def __init__( ) -> None: self._initial_soc = float(initial_soc) - try: - self._v_idx = self._pybamm_output_vars.index(_TERMINAL_VOLTAGE_VAR) - except ValueError: + if not any(v in self._pybamm_output_vars for v in _VOLTAGE_VAR_CANDIDATES): raise TypeError( - f"{type(self).__name__}._pybamm_output_vars must contain " - f"'{_TERMINAL_VOLTAGE_VAR}'." - ) from None + f"{type(self).__name__}._pybamm_output_vars must contain one of " + f"{_VOLTAGE_VAR_CANDIDATES}." + ) if model is None: model = pybamm.lithium_ion.SPMe( options={"thermal": self._thermal_option, **self._thermal_extra_options} ) + else: + model = _inject_thermal_options( + model, + {"thermal": self._thermal_option, **self._thermal_extra_options}, + ) self._parameter_values = _prepare_parameter_values(parameter_values) try: @@ -113,9 +287,35 @@ def __init__( parameter_values=self._parameter_values, solver=pybamm_solver, ) - sim.build(initial_soc=self._initial_soc, inputs=_DEFAULT_INPUTS) + _build_simulation(sim, model, self._initial_soc) + + available = sim.built_model.variables - all_out_vars = self._pybamm_output_vars + ["Discharge capacity [A.h]"] + # 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 + ) + self._v_idx = resolved_output_vars.index( + next(v for v in _VOLTAGE_VAR_CANDIDATES if v in resolved_output_vars) + ) + + extra_soc_var = soc_cap_var if soc_cap_var is not None else soc_direct_var + all_out_vars = resolved_output_vars + [extra_soc_var] casadi_objs = sim.built_model.export_casadi_objects( all_out_vars, input_parameter_order=list(_DEFAULT_INPUTS.keys()), @@ -123,17 +323,8 @@ def __init__( t_sym = casadi_objs["t"] x_sym = casadi_objs["x"] - z_sym = casadi_objs["z"] p_sym = casadi_objs["inputs"] - if z_sym.numel() > 0: - raise NotImplementedError( - f"{type(self).__name__}: the supplied PyBaMM model has " - f"{z_sym.numel()} algebraic variable(s) after discretisation " - "(DAE system). Only pure ODE models are supported. " - "Use SPMe or SPM instead of DFN." - ) - 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"]] @@ -153,7 +344,7 @@ def __init__( q_nominal = self._q_nominal initial_soc_val = float(initial_soc) - pybamm_output_vars = list(self._pybamm_output_vars) + soc_direct_scale = _detect_soc_direct_scale(sim, soc_direct_var) def _pack(u): return casadi.DM([float(u[0]), float(u[1])]) @@ -175,9 +366,13 @@ def jac_dyn(x, u, t): 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 pybamm_output_vars] - q_dis = float(out_var_fns["Discharge capacity [A.h]"](t, xv, p)) - soc = max(0.0, min(1.0, initial_soc_val - q_dis / q_nominal)) + 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: @@ -214,8 +409,10 @@ class _CoSimCellBase(Wrapper): differential-algebraic solve internally. Subclasses set ``_thermal_option``, ``_pybamm_output_vars`` and port labels. - ``_pybamm_output_vars`` must contain ``_TERMINAL_VOLTAGE_VAR``; its position - in the list is found dynamically. + ``_pybamm_output_vars`` must contain a terminal-voltage entry — either the + canonical ``"Terminal voltage [V]"`` or any alias in + ``_VOLTAGE_VAR_CANDIDATES``; the actual name is resolved against the built + model at construction time. """ _thermal_option: str = "" @@ -232,13 +429,12 @@ def __init__( ) -> None: self._initial_soc = float(initial_soc) - try: - self._v_idx = self._pybamm_output_vars.index(_TERMINAL_VOLTAGE_VAR) - except ValueError: + if not any(v in self._pybamm_output_vars for v in _VOLTAGE_VAR_CANDIDATES): raise TypeError( - f"{type(self).__name__}._pybamm_output_vars must contain " - f"'{_TERMINAL_VOLTAGE_VAR}'." - ) from None + f"{type(self).__name__}._pybamm_output_vars must contain one of " + f"{_VOLTAGE_VAR_CANDIDATES}." + ) + self._dt = float(dt) if self._dt <= 0.0: raise ValueError("dt must be positive") @@ -247,6 +443,11 @@ def __init__( model = pybamm.lithium_ion.SPMe( options={"thermal": self._thermal_option, **self._thermal_extra_options} ) + else: + model = _inject_thermal_options( + model, + {"thermal": self._thermal_option, **self._thermal_extra_options}, + ) self._model = model self._parameter_values = _prepare_parameter_values(parameter_values) @@ -264,6 +465,19 @@ def __init__( self._sim = self._build_sim() + available = self._sim.built_model.variables + self._resolved_output_vars, self._soc_cap_var, self._soc_direct_var = ( + _resolve_output_vars(self._pybamm_output_vars, available) + ) + self._v_idx = self._resolved_output_vars.index( + next(v for v in _VOLTAGE_VAR_CANDIDATES if v in self._resolved_output_vars) + ) + # Scale factor: 1.0 if SOC direct variable is in [0,1] fraction form, + # 1/100 if it is in percentage form (e.g. lead_acid "State of Charge"). + self._soc_direct_scale = _detect_soc_direct_scale( + self._sim, self._soc_direct_var + ) + self._last_outputs: npt.NDArray[np.float64] = self._initial_outputs() super().__init__(func=self._discrete_step, T=self._dt, tau=self._dt) @@ -278,7 +492,7 @@ def _build_sim(self) -> pybamm.Simulation: parameter_values=self._parameter_values, solver=self._pybamm_solver, ) - sim.build(initial_soc=self._initial_soc, inputs=_DEFAULT_INPUTS) + _build_simulation(sim, self._model, self._initial_soc) return sim def _initial_outputs(self) -> npt.NDArray[np.float64]: @@ -292,7 +506,10 @@ def _initial_outputs(self) -> npt.NDArray[np.float64]: not account for a non-default initial temperature or a non-zero current at t=0. """ - all_out_vars = self._pybamm_output_vars + ["Discharge capacity [A.h]"] + extra_soc_var = ( + self._soc_cap_var if self._soc_cap_var is not None else self._soc_direct_var + ) + all_out_vars = self._resolved_output_vars + [extra_soc_var] casadi_objs = self._sim.built_model.export_casadi_objects( all_out_vars, input_parameter_order=list(_DEFAULT_INPUTS.keys()), @@ -308,19 +525,22 @@ def _initial_outputs(self) -> npt.NDArray[np.float64]: z0 = casadi.Function("z0", [p_sym], [casadi_objs["z0"]])(p0) outputs: list[float] = [] - for name in self._pybamm_output_vars: + for name in self._resolved_output_vars: fn = casadi.Function( "v", [t_sym, x_sym, z_sym, p_sym], [casadi_objs["variables"][name]] ) outputs.append(float(fn(0.0, x0, z0, p0))) - q_dis_fn = casadi.Function( - "q", + extra_fn = casadi.Function( + "e", [t_sym, x_sym, z_sym, p_sym], - [casadi_objs["variables"]["Discharge capacity [A.h]"]], + [casadi_objs["variables"][extra_soc_var]], ) - q_dis = float(q_dis_fn(0.0, x0, z0, p0)) - soc = max(0.0, min(1.0, self._initial_soc - q_dis / self._q_nominal)) + extra_val = float(extra_fn(0.0, x0, z0, p0)) + if self._soc_cap_var is not None: + soc = max(0.0, min(1.0, self._initial_soc - extra_val / self._q_nominal)) + else: + soc = max(0.0, min(1.0, extra_val * self._soc_direct_scale)) outputs.append(soc) return np.array(outputs, dtype=np.float64) @@ -333,9 +553,13 @@ def _discrete_step(self, current: float, t_amb: float) -> npt.NDArray[np.float64 self._sim.step(dt=self._dt, inputs=inputs, save=False) sol = self._sim.solution - outputs = [float(sol[n].entries[-1]) for n in self._pybamm_output_vars] - q_dis = float(sol["Discharge capacity [A.h]"].entries[-1]) - soc = max(0.0, min(1.0, self._initial_soc - q_dis / self._q_nominal)) + outputs = [float(sol[n].entries[-1]) for n in self._resolved_output_vars] + if self._soc_cap_var is not None: + q_dis = float(sol[self._soc_cap_var].entries[-1]) + soc = max(0.0, min(1.0, self._initial_soc - q_dis / self._q_nominal)) + else: + raw = float(sol[self._soc_direct_var].entries[-1]) + soc = max(0.0, min(1.0, raw * self._soc_direct_scale)) outputs.append(soc) self._last_outputs = np.array(outputs, dtype=np.float64) diff --git a/tests/cells/_helpers.py b/tests/cells/_helpers.py new file mode 100644 index 0000000..a4c4a7e --- /dev/null +++ b/tests/cells/_helpers.py @@ -0,0 +1,201 @@ +"""Shared runner helpers and assertion helpers for cell test modules. + +Import these plain functions directly in test files; they are not pytest +fixtures. Every runner accepts an optional ``pybamm_solver`` keyword so +individual tests can override the PyBaMM solver (e.g. for ODE models that +require CasadiSolver instead of the default IDAKLUSolver). +""" + +from pathsim import Connection, Simulation +from pathsim.blocks import Constant +from pathsim.solvers import ESDIRK43 + +from pathsim_batt.cells import ( + CellCoSimElectrical, + CellCoSimElectrothermal, + CellElectrical, + CellElectrothermal, +) + +# --------------------------------------------------------------------------- +# Runner helpers +# --------------------------------------------------------------------------- + + +def run_electrical( + model, + pv, + current=1.0, + t_cell=298.15, + duration=2, + dt=1.0, + initial_soc=1.0, + pybamm_solver=None, +): + kwargs = {} + if pybamm_solver is not None: + kwargs["pybamm_solver"] = pybamm_solver + cell = CellElectrical( + model=model, parameter_values=pv, initial_soc=initial_soc, **kwargs + ) + I_src = Constant(current) + T_src = Constant(t_cell) + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell["T_cell"]), + ], + dt=dt, + Solver=ESDIRK43, + ) + sim.run(duration) + return cell + + +def run_electrothermal( + model, + pv, + current=1.0, + t_amb=298.15, + duration=2, + dt=1.0, + initial_soc=1.0, + pybamm_solver=None, +): + kwargs = {} + if pybamm_solver is not None: + kwargs["pybamm_solver"] = pybamm_solver + cell = CellElectrothermal( + model=model, parameter_values=pv, initial_soc=initial_soc, **kwargs + ) + I_src = Constant(current) + T_src = Constant(t_amb) + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell["T_amb"]), + ], + dt=dt, + Solver=ESDIRK43, + ) + sim.run(duration) + return cell + + +def run_cosim_electrical( + model, + pv, + current=1.0, + t_cell=298.15, + duration=2, + cosim_dt=1.0, + initial_soc=1.0, + pybamm_solver=None, +): + kwargs = {} + if pybamm_solver is not None: + kwargs["pybamm_solver"] = pybamm_solver + cell = CellCoSimElectrical( + model=model, + parameter_values=pv, + dt=cosim_dt, + initial_soc=initial_soc, + **kwargs, + ) + I_src = Constant(current) + T_src = Constant(t_cell) + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell["T_cell"]), + ], + dt=cosim_dt / 2, + Solver=ESDIRK43, + ) + sim.run(duration) + return cell + + +def run_cosim_electrothermal( + model, + pv, + current=1.0, + t_amb=298.15, + duration=2, + cosim_dt=1.0, + initial_soc=1.0, + pybamm_solver=None, +): + kwargs = {} + if pybamm_solver is not None: + kwargs["pybamm_solver"] = pybamm_solver + cell = CellCoSimElectrothermal( + model=model, + parameter_values=pv, + dt=cosim_dt, + initial_soc=initial_soc, + **kwargs, + ) + I_src = Constant(current) + T_src = Constant(t_amb) + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell["T_amb"]), + ], + dt=cosim_dt / 2, + Solver=ESDIRK43, + ) + sim.run(duration) + return cell + + +# --------------------------------------------------------------------------- +# Assertion helpers +# --------------------------------------------------------------------------- + + +def assert_electrical_outputs(test, cell, v_lo, v_hi, check_q_dot_nonneg=True): + """Assert physically plausible outputs for a CellElectrical run. + + Checks: + - V in (v_lo, v_hi) + - Q_dot >= 0 (optional) + - SOC in (0, 1] + """ + V = float(cell.outputs[0]) + Q = float(cell.outputs[1]) + soc = float(cell.outputs[2]) + test.assertGreater(V, v_lo, f"V={V:.3f} below lower cutoff {v_lo}") + test.assertLess(V, v_hi, f"V={V:.3f} above upper cutoff {v_hi}") + if check_q_dot_nonneg: + test.assertGreaterEqual(Q, 0.0, f"Q_dot={Q:.4f} is negative") + test.assertGreater(soc, 0.0) + test.assertLessEqual(soc, 1.0) + + +def assert_electrothermal_outputs(test, cell, v_lo, v_hi, check_q_dot_nonneg=True): + """Assert physically plausible outputs for a CellElectrothermal run. + + Checks: + - V in (v_lo, v_hi) + - T in (250, 400) K + - Q_dot >= 0 (optional) + - SOC in (0, 1] + """ + V = float(cell.outputs[0]) + T = float(cell.outputs[1]) + Q = float(cell.outputs[2]) + soc = float(cell.outputs[3]) + test.assertGreater(V, v_lo, f"V={V:.3f} below lower cutoff {v_lo}") + test.assertLess(V, v_hi, f"V={V:.3f} above upper cutoff {v_hi}") + test.assertGreater(T, 250.0, f"T={T:.1f} K unreasonably cold") + test.assertLess(T, 400.0, f"T={T:.1f} K unreasonably hot") + if check_q_dot_nonneg: + test.assertGreaterEqual(Q, 0.0, f"Q_dot={Q:.4f} is negative") + test.assertGreater(soc, 0.0) + test.assertLessEqual(soc, 1.0) diff --git a/tests/cells/test_ecm.py b/tests/cells/test_ecm.py new file mode 100644 index 0000000..985c0e6 --- /dev/null +++ b/tests/cells/test_ecm.py @@ -0,0 +1,168 @@ +"""Tests for equivalent_circuit.Thevenin (ECM) model. + +ECM_Example cutoffs: lower 3.2 V, upper 4.2 V, nominal capacity 100 A·h. + +ECM's reversible heat term can be negative (entropy effect), so Q_dot is +not constrained to be non-negative in the smoke tests. + +Co-simulation blocks are started at SOC=0.99 to avoid PyBaMM's internal +"Maximum SoC" event: the event value is zero at exactly SOC=1.0, which +the solver requires to be strictly positive at t=0. +""" + +import unittest + +import pybamm + +from pathsim_batt.cells import CellElectrical + +from ._helpers import ( + assert_electrical_outputs, + assert_electrothermal_outputs, + run_cosim_electrical, + run_cosim_electrothermal, + run_electrical, + run_electrothermal, +) + + +class TestECM(unittest.TestCase): + """equivalent_circuit.Thevenin with ECM_Example parameters (ODE — all blocks). + + ECM_Example cutoffs: lower 3.2 V, upper 4.2 V, nominal capacity 100 A·h. + + ECM's reversible heat term can be negative (entropy effect), so Q_dot is + not constrained to be non-negative in these tests. + + Co-simulation blocks are started at SOC=0.9 to avoid PyBaMM's internal + "Maximum SoC" event firing immediately at the upper boundary. + """ + + def setUp(self): + self.pv = pybamm.ParameterValues("ECM_Example") + self.v_lo = float(self.pv["Lower voltage cut-off [V]"]) + self.v_hi = float(self.pv["Upper voltage cut-off [V]"]) + + def _model(self): + return pybamm.equivalent_circuit.Thevenin() + + def test_electrical_smoke(self): + cell = run_electrical(self._model(), self.pv, current=10.0) + assert_electrical_outputs( + self, cell, self.v_lo, self.v_hi, check_q_dot_nonneg=False + ) + + def test_electrothermal_smoke(self): + cell = run_electrothermal(self._model(), self.pv, current=10.0) + assert_electrothermal_outputs( + self, cell, self.v_lo, self.v_hi, check_q_dot_nonneg=False + ) + + def test_cosim_electrical_smoke(self): + # Start at 0.99: PyBaMM's Maximum-SoC event fires at exactly SOC=1.0. + cell = run_cosim_electrical( + self._model(), self.pv, current=10.0, initial_soc=0.99 + ) + assert_electrical_outputs( + self, cell, self.v_lo, self.v_hi, check_q_dot_nonneg=False + ) + + def test_cosim_electrothermal_smoke(self): + cell = run_cosim_electrothermal( + self._model(), self.pv, current=10.0, initial_soc=0.99 + ) + assert_electrothermal_outputs( + self, cell, self.v_lo, self.v_hi, check_q_dot_nonneg=False + ) + + def test_electrical_soc_decreases(self): + cell = run_electrical(self._model(), self.pv, current=10.0, duration=60) + self.assertLess(float(cell.outputs[2]), 1.0) + + 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) + + def test_initial_soc_reflected_in_output(self): + """Output SOC must match initial_soc after zero-current run of 1 s. + + Verifies that set_initial_state is correctly wired to the CasADi output. + """ + cell = run_electrical( + self._model(), + self.pv, + current=0.0, + initial_soc=0.8, + duration=1, + ) + self.assertAlmostEqual( + float(cell.outputs[2]), + 0.8, + delta=0.01, + msg="Output SOC does not reflect initial_soc=0.8", + ) + + def test_soc_decrease_magnitude(self): + """Actual ΔSOC must match Coulombic prediction within 5 %. + + ECM's SoC state is purely Coulombic, so this should be tight. + """ + current = 10.0 + duration = 360 + pv = self.pv + q_nominal = float(pv["Nominal cell capacity [A.h]"]) + expected_delta = current * duration / 3600.0 / q_nominal + cell = run_electrical( + self._model(), + pv, + current=current, + initial_soc=1.0, + duration=duration, + ) + actual_delta = 1.0 - float(cell.outputs[2]) + self.assertAlmostEqual( + actual_delta, + expected_delta, + delta=expected_delta * 0.05, + msg=( + f"ΔSOC={actual_delta:.5f} deviates from Coulombic " + f"prediction {expected_delta:.5f} by more than 5 %" + ), + ) + + def test_q_dot_nonzero_during_discharge(self): + """abs(Q_dot) must be non-zero while current flows. + + ECM Q_dot can be negative (reversible heat) but must not be zero + during current flow. + """ + cell = run_electrical(self._model(), self.pv, current=10.0, duration=10) + self.assertGreater( + abs(float(cell.outputs[1])), + 1e-6, + "Q_dot is zero during discharge — heat generation not wired", + ) + + def test_tamb_affects_temperature(self): + """A warmer ambient temperature must yield a higher output cell temperature.""" + cell_cold = run_electrothermal( + self._model(), self.pv, current=10.0, t_amb=278.15, duration=120 + ) + cell_warm = run_electrothermal( + self._model(), self.pv, current=10.0, t_amb=318.15, duration=120 + ) + T_cold = float(cell_cold.outputs[1]) + T_warm = float(cell_warm.outputs[1]) + self.assertLess( + T_cold, + T_warm, + msg=( + f"Warmer ambient did not yield higher cell temperature: " + f"T_cold={T_cold:.2f} K, T_warm={T_warm:.2f} K" + ), + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/cells/test_lead_acid.py b/tests/cells/test_lead_acid.py new file mode 100644 index 0000000..0d06184 --- /dev/null +++ b/tests/cells/test_lead_acid.py @@ -0,0 +1,208 @@ +"""Tests for lead_acid model families: LOQS (ODE) and Full (DAE). + +Block / model matrix covered +----------------------------- +lead_acid.LOQS — ODE → all 4 blocks +lead_acid.Full — DAE → CoSim blocks only +""" + +import unittest + +import pybamm +from pathsim import Connection, Simulation +from pathsim.blocks import Constant +from pathsim.solvers import ESDIRK43 + +from pathsim_batt.cells import ( + CellCoSimElectrical, + CellCoSimElectrothermal, + CellElectrical, + CellElectrothermal, +) + +from ._helpers import ( + assert_electrical_outputs, + assert_electrothermal_outputs, + run_cosim_electrical, + run_cosim_electrothermal, + run_electrical, + run_electrothermal, +) + +# --------------------------------------------------------------------------- +# lead_acid.LOQS (ODE — all 4 blocks) +# --------------------------------------------------------------------------- + + +class TestLeadAcidLOQS(unittest.TestCase): + """lead_acid.LOQS with Sulzer2019 parameters (ODE model — all blocks). + + Sulzer2019 cutoffs: lower 1.75 V, upper 2.42 V, nominal capacity 17 A·h. + """ + + def setUp(self): + self.pv = pybamm.ParameterValues("Sulzer2019") + self.v_lo = float(self.pv["Lower voltage cut-off [V]"]) + self.v_hi = float(self.pv["Upper voltage cut-off [V]"]) + + def _model(self): + return pybamm.lead_acid.LOQS() + + 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) + + 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_cosim_electrical_smoke(self): + # LOQS is an ODE; IDAKLUSolver (co-sim default) requires a Jacobian + # for ODE models and errors, so use CasadiSolver explicitly. + solver = pybamm.CasadiSolver(mode="safe") + cell = CellCoSimElectrical( + model=self._model(), + parameter_values=self.pv, + pybamm_solver=solver, + dt=1.0, + ) + I_src = Constant(17.0) + T_src = Constant(298.15) + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell["T_cell"]), + ], + dt=0.5, + Solver=ESDIRK43, + ) + sim.run(2) + assert_electrical_outputs(self, cell, self.v_lo, self.v_hi) + + def test_cosim_electrothermal_smoke(self): + solver = pybamm.CasadiSolver(mode="safe") + cell = CellCoSimElectrothermal( + model=self._model(), + parameter_values=self.pv, + pybamm_solver=solver, + dt=1.0, + ) + I_src = Constant(17.0) + T_src = Constant(298.15) + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell["T_amb"]), + ], + dt=0.5, + Solver=ESDIRK43, + ) + sim.run(2) + assert_electrothermal_outputs(self, cell, self.v_lo, self.v_hi) + + 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) + + 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) + + def test_q_dot_nonzero_during_discharge(self): + """Q_dot must be strictly positive during discharge (isothermal LOQS). + + The block automatically injects the heat-source calculation flag into + isothermal models that lack it, so no manual option is needed. + """ + cell = run_electrical(self._model(), self.pv, current=17.0, duration=60) + self.assertGreater( + float(cell.outputs[1]), + 0.0, + "Q_dot is zero — thermal model may not compute heat sources", + ) + + def test_tamb_affects_temperature(self): + """A warmer ambient temperature must yield a higher output cell temperature.""" + solver = pybamm.CasadiSolver(mode="safe") + cell_cold = run_electrothermal( + self._model(), + self.pv, + current=17.0, + t_amb=278.15, + duration=120, + pybamm_solver=solver, + ) + cell_warm = run_electrothermal( + self._model(), + self.pv, + current=17.0, + t_amb=318.15, + duration=120, + pybamm_solver=solver, + ) + T_cold = float(cell_cold.outputs[1]) + T_warm = float(cell_warm.outputs[1]) + self.assertLess( + T_cold, + T_warm, + msg=( + f"Warmer ambient (318.15 K) did not yield higher cell temperature: " + f"T_cold={T_cold:.2f} K, T_warm={T_warm:.2f} K" + ), + ) + + def test_soc_scale_factor(self): + """SOC must be well below 1.0 after sustained discharge. + + Guards against PyBaMM's percentage-form 'State of Charge' being + misidentified as a fraction, which would clamp all values to 1.0. + """ + cell = run_electrical(self._model(), self.pv, current=17.0, duration=360) + self.assertLess( + float(cell.outputs[2]), + 0.95, + "SOC did not decrease — possible percentage/fraction scale error", + ) + + +# --------------------------------------------------------------------------- +# lead_acid.Full (DAE — co-simulation only) +# --------------------------------------------------------------------------- + + +class TestLeadAcidFull(unittest.TestCase): + """lead_acid.Full with Sulzer2019 parameters (DAE model — co-sim only).""" + + def setUp(self): + self.pv = pybamm.ParameterValues("Sulzer2019") + self.v_lo = float(self.pv["Lower voltage cut-off [V]"]) + self.v_hi = float(self.pv["Upper voltage cut-off [V]"]) + + 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_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_cosim_electrical_smoke(self): + cell = run_cosim_electrical(self._model(), self.pv, current=17.0) + assert_electrical_outputs(self, cell, self.v_lo, self.v_hi) + + def test_cosim_electrothermal_smoke(self): + cell = run_cosim_electrothermal(self._model(), self.pv, current=17.0) + assert_electrothermal_outputs(self, cell, self.v_lo, self.v_hi) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/cells/test_pybamm_cell.py b/tests/cells/test_pybamm_cell.py index 321e7d2..0db97cf 100644 --- a/tests/cells/test_pybamm_cell.py +++ b/tests/cells/test_pybamm_cell.py @@ -248,6 +248,51 @@ def _run_and_get_voltage(T_cell): msg="T_cell input has no effect on terminal voltage", ) + def test_initial_soc_reflected_in_output(self): + """Output SOC must match initial_soc=0.5 after 1 s at zero current.""" + cell = CellElectrical(initial_soc=0.5) + sim = self._make_simulation(cell, current=0.0, T_cell=298.15) + sim.run(1) + self.assertAlmostEqual( + float(cell.outputs[2]), + 0.5, + delta=0.01, + msg="Output SOC does not reflect initial_soc=0.5", + ) + + def test_soc_decrease_magnitude(self): + """Actual ΔSOC must match Coulombic prediction within 2 %. + + SOC is computed from discharge capacity which is exactly Coulombic. + """ + current = 5.0 + duration = 360 + dt = 10.0 + cell = CellElectrical(initial_soc=1.0) + I_src = Constant(current) + T_src = Constant(298.15) + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell["T_cell"]), + ], + dt=dt, + Solver=ESDIRK43, + ) + sim.run(duration) + expected_delta = current * duration / 3600.0 / cell._q_nominal + actual_delta = 1.0 - float(cell.outputs[2]) + self.assertAlmostEqual( + actual_delta, + expected_delta, + delta=expected_delta * 0.02, + msg=( + f"ΔSOC={actual_delta:.5f} deviates from Coulombic " + f"prediction {expected_delta:.5f} by more than 2 %" + ), + ) + class TestElectrothermal(unittest.TestCase): """Integration tests for CellElectrothermal — PathSim integrates the PyBaMM ODE.""" @@ -452,6 +497,23 @@ def _run_and_get_voltage(T_cell): msg="T_cell input has no effect on terminal voltage", ) + def test_initial_soc_reflected_in_output(self): + """Output SOC must match initial_soc=0.5 after 1 s (shorter than one step). + + When the elapsed time is shorter than the co-sim macro-step dt, the + block has not yet fired its discrete step and outputs remain at the + initial values — so SOC should equal the configured initial_soc. + """ + cell = CellCoSimElectrical(initial_soc=0.5, dt=10.0) + sim = self._make_simulation(cell, current=0.0, T_cell=298.15) + sim.run(1) + self.assertAlmostEqual( + float(cell.outputs[2]), + 0.5, + delta=0.01, + msg="Output SOC does not reflect initial_soc=0.5", + ) + class TestCoSimulationElectrothermal(unittest.TestCase): """Integration tests for CellCoSimElectrothermal — PyBaMM performs the stepping."""