diff --git a/src/pathsim_batt/cells/pybamm_cell.py b/src/pathsim_batt/cells/pybamm_cell.py index 40ce974..551633d 100644 --- a/src/pathsim_batt/cells/pybamm_cell.py +++ b/src/pathsim_batt/cells/pybamm_cell.py @@ -14,6 +14,7 @@ import numpy.typing as npt import pybamm from pathsim.blocks import DynamicalSystem, Wrapper +from pathsim.events import ZeroCrossingDown # HELPERS ============================================================================= @@ -80,6 +81,12 @@ def __init__( ) self._parameter_values = _prepare_parameter_values(parameter_values) + self._v_lower = float( + self._parameter_values["Lower voltage cut-off [V]"] + ) + self._v_upper = float( + self._parameter_values["Upper voltage cut-off [V]"] + ) pybamm_solver = pybamm_solver or pybamm.CasadiSolver(mode="safe") @@ -166,6 +173,57 @@ def func_alg(x, u, t): def __len__(self) -> int: return len(self._pybamm_output_vars) + 1 + def termination_events(self, sim) -> list: + """Return PathSim events that mirror PyBaMM's voltage cut-off conditions. + + Registers two :class:`~pathsim.events.ZeroCrossingDown` events on this + cell's terminal voltage output (port ``V``, index 0): + + * **Under-voltage** — fires when ``V`` falls to ``Lower voltage + cut-off [V]`` (identical to PyBaMM's ``Minimum voltage [V]`` + termination event). + * **Over-voltage** — fires when ``V`` rises to ``Upper voltage + cut-off [V]`` (identical to PyBaMM's ``Maximum voltage [V]`` + termination event). + + Each event calls ``sim.stop()`` when it resolves, halting the PathSim + time-stepping loop at the crossing point. + + Parameters + ---------- + sim : pathsim.Simulation + The simulation this cell is part of. The events call + ``sim.stop()`` on resolution. + + Returns + ------- + list[ZeroCrossingDown] + Two events; add them to the simulation with + ``sim.add_event(e)`` or ``for e in cell.termination_events(sim): sim.add_event(e)``. + + Example + ------- + .. code-block:: python + + sim = Simulation(blocks=[...], connections=[...], dt=1.0, Solver=ESDIRK43) + for event in cell.termination_events(sim): + sim.add_event(event) + sim.run(3600) + """ + v_lower = self._v_lower + v_upper = self._v_upper + outputs = self.outputs + + under_voltage = ZeroCrossingDown( + func_evt=lambda t: float(outputs[0]) - v_lower, + func_act=lambda t: sim.stop(), + ) + over_voltage = ZeroCrossingDown( + func_evt=lambda t: v_upper - float(outputs[0]), + func_act=lambda t: sim.stop(), + ) + return [under_voltage, over_voltage] + def reset(self) -> None: super().reset() @@ -206,6 +264,12 @@ def __init__( self._model = model self._parameter_values = _prepare_parameter_values(parameter_values) + self._v_lower = float( + self._parameter_values["Lower voltage cut-off [V]"] + ) + self._v_upper = float( + self._parameter_values["Upper voltage cut-off [V]"] + ) self._pybamm_solver = pybamm_solver or pybamm.IDAKLUSolver() self._q_nominal = float(self._parameter_values["Nominal cell capacity [A.h]"]) @@ -213,6 +277,11 @@ def __init__( self._last_outputs: npt.NDArray[np.float64] = self._initial_outputs() + # registered stop callbacks — populated by termination_events() + self._term_callbacks: list = [] + # flag set by _discrete_step so update() knows when fresh outputs exist + self._just_stepped: bool = False + super().__init__(func=self._discrete_step, T=self._dt, tau=self._dt) # ensure outputs are valid before first scheduled sample @@ -229,15 +298,39 @@ def _build_sim(self) -> pybamm.Simulation: return sim def _initial_outputs(self) -> npt.NDArray[np.float64]: - """Return placeholder outputs for t=0 before the first solver step. - - The co-simulation takes its first real sample at t=dt, so this - placeholder is only held for one macro-step. All outputs are zero - except SOC, which is set to the user-supplied initial value. + """Compute physically correct outputs at t=0 from the built PyBaMM model. + + Uses the same CasADi export approach as ``_CellBase`` to evaluate each + output variable at the initial state vector and zero current, so that + the terminal voltage placeholder is the true open-circuit voltage. This + ensures that voltage-threshold events (``termination_events()``) have a + correct, positive starting value and can detect the first downward + crossing correctly. """ - out = np.zeros(len(self._pybamm_output_vars) + 1, dtype=np.float64) - out[-1] = self._initial_soc # SOC is always the last output - return out + all_out_vars = self._pybamm_output_vars + ["Discharge capacity [A.h]"] + casadi_objs = self._sim.built_model.export_casadi_objects( + all_out_vars, + input_parameter_order=list(_DEFAULT_INPUTS.keys()), + ) + t_sym = casadi_objs["t"] + x_sym = casadi_objs["x"] + p_sym = casadi_objs["inputs"] + p0 = casadi.DM(list(_DEFAULT_INPUTS.values())) + x0 = casadi.Function("x0", [p_sym], [casadi_objs["x0"]])(p0) + + outputs: list[float] = [] + for name in self._pybamm_output_vars: + fn = casadi.Function("v", [t_sym, x_sym, p_sym], [casadi_objs["variables"][name]]) + outputs.append(float(fn(0.0, x0, p0))) + + q_dis_fn = casadi.Function( + "q", [t_sym, x_sym, p_sym], [casadi_objs["variables"]["Discharge capacity [A.h]"]] + ) + q_dis = float(q_dis_fn(0.0, x0, p0)) + soc = max(0.0, min(1.0, self._initial_soc - q_dis / self._q_nominal)) + outputs.append(soc) + + return np.array(outputs, dtype=np.float64) def _discrete_step(self, current: float, t_amb: float) -> npt.NDArray[np.float64]: inputs = { @@ -253,18 +346,82 @@ def _discrete_step(self, current: float, t_amb: float) -> npt.NDArray[np.float64 outputs.append(soc) self._last_outputs = np.array(outputs, dtype=np.float64) + self._just_stepped = True return self._last_outputs def update(self, t: float) -> None: - self.outputs.update_from_array(self._last_outputs) + """Check voltage cut-off callbacks after each PyBaMM macro-step. + + PathSim calls ``update()`` on every block after each event resolves + (including after the internal :class:`~pathsim.events.Schedule` event + fires and updates outputs). The ``_just_stepped`` flag distinguishes + this post-step call from the earlier pre-event call where outputs are + still stale, ensuring the termination check only runs once per + macro-step and only after fresh outputs are available. + """ + if self._just_stepped: + self._just_stepped = False + V = float(self.outputs[0]) + for cb in self._term_callbacks: + cb(V) def __len__(self) -> int: return len(self._pybamm_output_vars) + 1 + def termination_events(self, sim) -> list: + """Return PathSim events that mirror PyBaMM's voltage cut-off conditions. + + For co-simulation blocks PyBaMM's own solver clamps the terminal voltage + at the cut-off and stops advancing internally, but never signals PathSim. + This method registers **post-step callbacks** (called from + :meth:`update` after each :class:`~pathsim.events.Schedule` macro-step + fires) that call ``sim.stop()`` as soon as the clamped output is + detected. Two :class:`~pathsim.events.ZeroCrossingDown` events are + also returned for API consistency with + :class:`_CellBase.termination_events` — add them to the simulation + with ``sim.add_event(e)`` as a belt-and-suspenders complement for + adaptive-stepping scenarios. + + Parameters + ---------- + sim : pathsim.Simulation + The simulation this cell is part of. + + Returns + ------- + list[ZeroCrossingDown] + Two events (under-voltage, over-voltage). + """ + v_lower = self._v_lower + v_upper = self._v_upper + + def _under_cb(V: float) -> None: + if V <= v_lower: + sim.stop() + + def _over_cb(V: float) -> None: + if V >= v_upper: + sim.stop() + + self._term_callbacks.extend([_under_cb, _over_cb]) + + # Belt-and-suspenders PathSim events (API parity with _CellBase). + outputs = self.outputs + under_voltage = ZeroCrossingDown( + func_evt=lambda t: float(outputs[0]) - v_lower, + func_act=lambda t: sim.stop(), + ) + over_voltage = ZeroCrossingDown( + func_evt=lambda t: v_upper - float(outputs[0]), + func_act=lambda t: sim.stop(), + ) + return [under_voltage, over_voltage] + def reset(self) -> None: super().reset() self._sim = self._build_sim() self._last_outputs = self._initial_outputs() + self._just_stepped = False self.outputs.update_from_array(self._last_outputs) diff --git a/tests/cells/test_pybamm_cell.py b/tests/cells/test_pybamm_cell.py index 681efe1..9a30182 100644 --- a/tests/cells/test_pybamm_cell.py +++ b/tests/cells/test_pybamm_cell.py @@ -555,5 +555,194 @@ def _run_and_get_T_cell(T_amb): ) +class TestTerminationEvents(unittest.TestCase): + """Tests that termination_events() stops the simulation at voltage cut-offs. + + Both the non-CoSim (``_CellBase``) and CoSim (``_CoSimCellBase``) families + are tested. A very low initial SOC combined with a high discharge current + ensures the lower voltage cut-off is reached quickly. Without the events + wired in, the non-CoSim cell would diverge to nonsensical negative voltages, + and the CoSim cell would continue forever at the clamped cut-off value. + + PyBaMM's ``Chen2020`` parameter set has: + * Lower voltage cut-off: 2.5 V + * Upper voltage cut-off: 4.2 V + """ + + # -- helpers --------------------------------------------------------------- + + def _run_with_events(self, cell, current, T_input_port, T_value, cosim_dt=None): + """Build a Simulation, wire termination events, and run for up to 600 s.""" + I_src = Constant(current) + T_src = Constant(T_value) + dt_ps = cosim_dt if cosim_dt is not None else 5.0 + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell[T_input_port]), + ], + dt=dt_ps, + Solver=ESDIRK43, + ) + for event in cell.termination_events(sim): + sim.add_event(event) + sim.run(600) + return sim + + def _run_without_events(self, cell, current, T_input_port, T_value, cosim_dt=None): + """Same but without registering termination events.""" + I_src = Constant(current) + T_src = Constant(T_value) + dt_ps = cosim_dt if cosim_dt is not None else 5.0 + sim = Simulation( + blocks=[I_src, T_src, cell], + connections=[ + Connection(I_src, cell["I"]), + Connection(T_src, cell[T_input_port]), + ], + dt=dt_ps, + Solver=ESDIRK43, + ) + sim.run(600) + return sim + + # -- non-CoSim (CellElectrical) -------------------------------------------- + + def test_non_cosim_stops_before_negative_voltage(self): + """With events registered, CellElectrical must stop before V < 0. + + Without events PathSim integrates through the physical cutoff and the + ODE produces nonsensical negative voltages. With events wired in the + simulation must halt at or just above the lower cut-off. + """ + cell = CellElectrical(initial_soc=0.02) + sim = self._run_with_events(cell, current=10.0, T_input_port="T_cell", T_value=298.15) + V_final = float(cell.outputs[0]) + # Simulation must have stopped early (well before 600 s) + self.assertLess(sim.time, 600.0, "simulation did not stop early") + # Voltage must be at or above the lower cut-off, never negative + self.assertGreaterEqual( + V_final, + cell._v_lower - 0.5, # small tolerance for step overshoot + f"terminal voltage {V_final:.3f} V is far below cut-off {cell._v_lower} V", + ) + self.assertGreater( + V_final, + 0.0, + f"terminal voltage went negative ({V_final:.3f} V) despite termination events", + ) + + def test_non_cosim_without_events_continues_past_cutoff(self): + """Without events, CellElectrical integrates past the cut-off (regression guard). + + This test confirms the *bug* that termination_events() is designed to fix: + PathSim does not stop on its own, and the voltage diverges wildly. + """ + cell = CellElectrical(initial_soc=0.02) + sim = self._run_without_events(cell, current=10.0, T_input_port="T_cell", T_value=298.15) + V_final = float(cell.outputs[0]) + # Without events the simulation runs the full 600 s + self.assertAlmostEqual(sim.time, 600.0, delta=10.0) + # And the voltage should have diverged to a physically impossible value + self.assertLess( + V_final, + 0.0, + "expected voltage divergence without termination events, got " + f"V={V_final:.3f} V — the bug may have been fixed elsewhere", + ) + + def test_non_cosim_termination_events_returns_two_events(self): + """termination_events() must return exactly two ZeroCrossingDown instances.""" + from pathsim.events import ZeroCrossingDown + + cell = CellElectrical() + I_src = Constant(1.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=1.0, + Solver=ESDIRK43, + ) + events = cell.termination_events(sim) + self.assertEqual(len(events), 2) + for e in events: + self.assertIsInstance(e, ZeroCrossingDown) + + def test_non_cosim_cutoff_values_match_parameter_values(self): + """_v_lower/_v_upper must match the Chen2020 parameter set.""" + pv = pybamm.ParameterValues("Chen2020") + cell = CellElectrical(parameter_values=pv) + self.assertAlmostEqual(cell._v_lower, float(pv["Lower voltage cut-off [V]"])) + self.assertAlmostEqual(cell._v_upper, float(pv["Upper voltage cut-off [V]"])) + + # -- CoSim (CellCoSimElectrical) ------------------------------------------- + + def test_cosim_stops_before_negative_voltage(self): + """With events registered, CellCoSimElectrical must stop at the cut-off. + + Without events PyBaMM freezes internally but PathSim keeps ticking + forever with the clamped output (2.5 V). With events wired in, + PathSim must stop as soon as the voltage first touches the cut-off. + """ + cell = CellCoSimElectrical(initial_soc=0.02, dt=10.0) + sim = self._run_with_events( + cell, current=10.0, T_input_port="T_cell", T_value=298.15, cosim_dt=10.0 + ) + V_final = float(cell.outputs[0]) + self.assertLess(sim.time, 600.0, "CoSim simulation did not stop early") + self.assertGreaterEqual( + V_final, + cell._v_lower - 0.5, + f"terminal voltage {V_final:.3f} V is far below cut-off {cell._v_lower} V", + ) + + def test_cosim_without_events_runs_full_duration(self): + """Without events, CellCoSimElectrical runs the full duration with frozen output. + + This is the CoSim-side symptom of the missing stop condition: PyBaMM + clamps at the cut-off and the simulation never terminates early. + """ + cell = CellCoSimElectrical(initial_soc=0.02, dt=10.0) + sim = self._run_without_events( + cell, current=10.0, T_input_port="T_cell", T_value=298.15, cosim_dt=10.0 + ) + # Simulation must have run the full 600 s (no early stop) + self.assertAlmostEqual(sim.time, 600.0, delta=10.0) + # Voltage must be frozen at the cut-off (PyBaMM clamps, does not diverge) + self.assertAlmostEqual( + float(cell.outputs[0]), + cell._v_lower, + delta=0.1, + msg="CoSim output should be clamped at lower cut-off without events", + ) + + def test_cosim_termination_events_returns_two_events(self): + """termination_events() must return exactly two ZeroCrossingDown instances.""" + from pathsim.events import ZeroCrossingDown + + cell = CellCoSimElectrical(dt=1.0) + I_src = Constant(1.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=1.0, + Solver=ESDIRK43, + ) + events = cell.termination_events(sim) + self.assertEqual(len(events), 2) + for e in events: + self.assertIsInstance(e, ZeroCrossingDown) + + def test_cosim_cutoff_values_match_parameter_values(self): + """_v_lower/_v_upper must match the Chen2020 parameter set (CoSim block).""" + pv = pybamm.ParameterValues("Chen2020") + cell = CellCoSimElectrical(parameter_values=pv, dt=1.0) + self.assertAlmostEqual(cell._v_lower, float(pv["Lower voltage cut-off [V]"])) + self.assertAlmostEqual(cell._v_upper, float(pv["Upper voltage cut-off [V]"])) + + if __name__ == "__main__": unittest.main()