From 9341e2f420b98d651a0bea501d2500a1dc3004ed Mon Sep 17 00:00:00 2001 From: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> Date: Mon, 19 Jan 2026 18:26:28 +0000 Subject: [PATCH 1/3] ParameterSubstitutor bugfixes (#5352) * fix: inputs dict bug w/ unprocessed child parameters, cache collision bugs * docs: changelog entry (cherry picked from commit 3517cff734cd3b35c3f6388a0ecf27e73af891bc) --- CHANGELOG.md | 2 ++ .../parameters/parameter_substitutor.py | 11 +------ .../test_parameters/test_parameter_values.py | 29 +++++++++++++++++++ 3 files changed, 32 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cf472a894a..3292bb548b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,8 @@ ## Bug fixes +- Fixed a bug where `ExpressionFunctionParameter` children were not properly substituted after JSON deserialization, causing `ElectrodeSOHSolver` to fail with "Variable not implemented" errors. ([#5352](https://github.com/pybamm-team/PyBaMM/pull/5352)) + # [v25.12.0](https://github.com/pybamm-team/PyBaMM/tree/v25.12.0) - 2026-01-15 ## Features diff --git a/src/pybamm/parameters/parameter_substitutor.py b/src/pybamm/parameters/parameter_substitutor.py index 65d435e625..bfc26af86f 100644 --- a/src/pybamm/parameters/parameter_substitutor.py +++ b/src/pybamm/parameters/parameter_substitutor.py @@ -297,7 +297,7 @@ def _process_function_parameter( inputs = { arg: child for arg, child in zip( - function_parameter.func_args, symbol.children, strict=False + function_parameter.func_args, new_children, strict=False ) } @@ -316,9 +316,6 @@ def _process_function_parameter( combined_store.update(inputs) combined_processor = ParameterSubstitutor(combined_store) - # Copy parent cache to combined processor - combined_processor._cache.update(self._cache) - # Process any FunctionParameter children first to avoid recursion for child in expression.pre_order(): if isinstance(child, pybamm.FunctionParameter): @@ -346,12 +343,6 @@ def _process_function_parameter( # Process function with combined processor to get a symbolic expression function = combined_processor.process_symbol(expression) - # Propagate newly cached symbols back to parent cache for reuse - # Only propagate symbols that don't depend on the local inputs - for sym, processed in combined_processor._cache.items(): - if sym not in self._cache and sym not in inputs.values(): - self._cache[sym] = processed - # Differentiate if necessary if symbol.diff_variable is None: # Use ones_like so that we get the right shapes diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index 2559c42ff1..d354f50ae6 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -1825,3 +1825,32 @@ def test_process_model_preserves_variables(self): assert var.id == original_var_ids[name], ( f"Variable '{name}' expression was modified by process_model" ) + + def test_json_roundtrip_with_electrode_soh_solver(self): + """Test json roundtrip with electrode soh_solver""" + # Serialize and deserialize parameter values + pv = pybamm.ParameterValues("Chen2020") + pv_dict = pv.to_json() + pv_reconstructed = pybamm.ParameterValues.from_json(pv_dict) + + # The ElectrodeSOHSolver uses OCP functions which are converted to + # ExpressionFunctionParameter after JSON roundtrip + param = pybamm.LithiumIonParameters() + esoh_solver = pybamm.lithium_ion.ElectrodeSOHSolver( + pv_reconstructed, direction=None, param=param + ) + + # Get inputs for the solver + Q_n = pv_reconstructed.evaluate(param.n.Q_init) + Q_p = pv_reconstructed.evaluate(param.p.Q_init) + Q_Li = pv_reconstructed.evaluate(param.Q_Li_particles_init) + inputs = {"Q_Li": Q_Li, "Q_n": Q_n, "Q_p": Q_p} + + # Should solve w/o issue + sol = esoh_solver.solve(inputs) + + # Verify the solution is valid + assert "x_100" in sol + assert "y_100" in sol + assert "Maximum theoretical energy [W.h]" in sol + assert sol["Maximum theoretical energy [W.h]"] > 0 From a02bdc4d8bbf2b2d6ad99094cb995642be4b0815 Mon Sep 17 00:00:00 2001 From: Robert Timms <43040151+rtimms@users.noreply.github.com> Date: Tue, 20 Jan 2026 14:56:32 +0000 Subject: [PATCH 2/3] Improve composite eSOH (#5347) * minor composite esoh fixes * style: pre-commit fixes * try split composite eSOH solve * use split for initial conditions on failure * style: pre-commit fixes * fix v_low/high vs soc_0/100 * consisten min/max sto calculations * update changelog * style: pre-commit fixes * PR comments * fix direction in solve_split * Update src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py Co-authored-by: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com> --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com> (cherry picked from commit 18461a186c0550722aa839f9bafd5c3d437e7b8c) --- CHANGELOG.md | 2 + .../lithium_ion/__init__.py | 6 +- .../lithium_ion/electrode_soh.py | 57 +- .../lithium_ion/electrode_soh_composite.py | 777 +++++++++++++++--- .../lithium_ion/electrode_soh_half_cell.py | 10 +- .../full_battery_models/lithium_ion/util.py | 67 +- .../test_lithium_ion/test_electrode_soh.py | 120 ++- 7 files changed, 899 insertions(+), 140 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3292bb548b..2637bbb731 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,8 @@ ## Features +- Improve the robustness of the `ElectrodeSOHComposite` model and ensure the calculation of equilibrium stoichiometries is consistent for models with hysteresis. ([#5347](https://github.com/pybamm-team/PyBaMM/pull/5347)) + ## Bug fixes - Fixed a bug where `ExpressionFunctionParameter` children were not properly substituted after JSON deserialization, causing `ElectrodeSOHSolver` to fail with "Variable not implemented" errors. ([#5352](https://github.com/pybamm-team/PyBaMM/pull/5352)) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py index 330b6c9e42..b8f106be7e 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py @@ -32,7 +32,11 @@ from .mpm import MPM from .msmr import MSMR from .basic_splitOCVR import SplitOCVR -from .util import check_if_composite +from .util import ( + check_if_composite, + get_equilibrium_direction, + get_lithiation_delithiation, +) __all__ = ['Yang2017', 'base_lithium_ion_model', 'basic_dfn', 'basic_dfn_composite', 'basic_dfn_half_cell', 'basic_spm', 'dfn', diff --git a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py index b769727d9e..982e22ab6f 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py @@ -7,7 +7,7 @@ import pybamm -from .util import _get_lithiation_delithiation +from .util import get_equilibrium_direction, get_lithiation_delithiation class _BaseElectrodeSOH(pybamm.BaseModel): @@ -96,6 +96,13 @@ class _ElectrodeSOH(_BaseElectrodeSOH): .. math:: y_0 = y_{100} + \\frac{Q}{Q_p}. + Stoichiometry limits (x_0, x_100, y_0, y_100) are evaluated in a + temperature-independent manner (at the reference temperature, ignoring entropy + effects), so that the SOC range is not temperature-dependent. In the presence of + a hysteresis model, equilibration in the charging direction is assumed for 100% + SOC (charging OCP branch), and equilibration in the discharging direction is + assumed for 0% SOC (discharging OCP branch). + """ def __init__( @@ -147,10 +154,22 @@ def __init__( x_100 = pybamm.InputParameter("x_100") y_100 = pybamm.InputParameter("y_100") Un_100 = Un( - x_100, T_ref, _get_lithiation_delithiation(direction, "negative", options) + x_100, + T_ref, + get_lithiation_delithiation( + get_equilibrium_direction("100", "negative", options), + "negative", + options, + ), ) Up_100 = Up( - y_100, T_ref, _get_lithiation_delithiation(direction, "positive", options) + y_100, + T_ref, + get_lithiation_delithiation( + get_equilibrium_direction("100", "positive", options), + "positive", + options, + ), ) # Define equations for 100% state of charge @@ -178,10 +197,22 @@ def __init__( var = y_100 y_0 = y_100 + Q / Q_p Un_0 = Un( - x_0, T_ref, _get_lithiation_delithiation(direction, "negative", options) + x_0, + T_ref, + get_lithiation_delithiation( + get_equilibrium_direction("0", "negative", options), + "negative", + options, + ), ) Up_0 = Up( - y_0, T_ref, _get_lithiation_delithiation(direction, "positive", options) + y_0, + T_ref, + get_lithiation_delithiation( + get_equilibrium_direction("0", "positive", options), + "positive", + options, + ), ) self.algebraic[var] = Up_0 - Un_0 - V_min self.initial_conditions[var] = pybamm.Scalar(0.1) @@ -658,14 +689,14 @@ def _check_esoh_feasible(self, inputs, direction): self.param.p.prim.U( y, T, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "positive", self.options ), ) - self.param.n.prim.U( x, T, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "negative", self.options ), ) @@ -784,14 +815,14 @@ def get_initial_stoichiometries( Up( y, T_init, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "positive", self.options ), ) - Un( x, T_init, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "negative", self.options ), ) @@ -902,7 +933,7 @@ def get_initial_ocps(self, initial_value, direction=None, tol=1e-6, inputs=None) self.param.n.prim.U( x, T_init, - _get_lithiation_delithiation(direction, "negative", self.options), + get_lithiation_delithiation(direction, "negative", self.options), ), inputs=inputs, ) @@ -910,7 +941,7 @@ def get_initial_ocps(self, initial_value, direction=None, tol=1e-6, inputs=None) self.param.p.prim.U( y, T_init, - _get_lithiation_delithiation(direction, "positive", self.options), + get_lithiation_delithiation(direction, "positive", self.options), ), inputs=inputs, ) @@ -963,12 +994,12 @@ def theoretical_energy_integral(self, inputs, points=1000, direction="discharge" self.param.p.prim.U( y_vals, T, - _get_lithiation_delithiation(direction, "positive", self.options), + get_lithiation_delithiation(direction, "positive", self.options), ) - self.param.n.prim.U( x_vals, T, - _get_lithiation_delithiation(direction, "negative", self.options), + get_lithiation_delithiation(direction, "negative", self.options), ), inputs=inputs, ).flatten() diff --git a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py index bb90205f85..1e18959f36 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_composite.py @@ -1,12 +1,80 @@ # # A model to calculate electrode-specific SOH, adapted for composite electrodes # +from __future__ import annotations + +from typing import Any + import pybamm -from .util import _get_lithiation_delithiation, check_if_composite +from .electrode_soh import _ElectrodeSOH +from .util import ( + check_if_composite, + get_equilibrium_direction, + get_lithiation_delithiation, +) + + +def _get_primary_only_options( + options: dict | pybamm.BatteryModelOptions | None, +) -> dict | None: + """ + Create options dict with only primary phase OCP settings. + + When composite options have tuple OCP settings like + (("single", "hysteresis"), "single"), the non-composite model needs only the + primary (first) element for each electrode to avoid incorrectly detecting + hysteresis from the secondary phase. + + Parameters + ---------- + options : dict or pybamm.BatteryModelOptions or None + Model options that may contain composite OCP settings. + + Returns + ------- + dict or None + Options dict with primary-only OCP settings, or None if input is None. + + Examples + -------- + >>> opts = _get_primary_only_options( + ... {"open-circuit potential": (("single", "current sigmoid"), "single")} + ... ) + >>> opts["open-circuit potential"] + ('single', 'single') + """ + if options is None: + return None + + if not isinstance(options, pybamm.BatteryModelOptions): + options = pybamm.BatteryModelOptions(options) + + options_dict = dict(options) + + ocp_option = options_dict.get("open-circuit potential") + if ocp_option is None: + return options_dict + + if isinstance(ocp_option, str): + return options_dict + + if isinstance(ocp_option, tuple) and len(ocp_option) == 2: + neg_ocp, pos_ocp = ocp_option + + if isinstance(neg_ocp, tuple): + neg_ocp = neg_ocp[0] + + if isinstance(pos_ocp, tuple): + pos_ocp = pos_ocp[0] + + options_dict["open-circuit potential"] = (neg_ocp, pos_ocp) + + return options_dict def _get_stoich_variables(options): + """Create stoichiometry variables for composite electrodes.""" variables = { "x_100_1": pybamm.Variable("x_100_1", bounds=(0, 1)), "y_100_1": pybamm.Variable("y_100_1", bounds=(0, 1)), @@ -29,6 +97,7 @@ def _get_stoich_variables(options): def _get_initial_conditions(options, soc_init): + """Get initial conditions for stoichiometry variables.""" variables = _get_stoich_variables(options) ics = {} # Avoid exact boundary values for better numerical stability @@ -50,6 +119,7 @@ def _get_initial_conditions(options, soc_init): def _get_direction(electrode): + """Get direction multiplier for electrode capacity calculations.""" if electrode == "positive": return pybamm.Scalar(-1) else: @@ -57,6 +127,7 @@ def _get_direction(electrode): def _get_prefix(electrode): + """Get stoichiometry variable prefix for electrode ('x' or 'y').""" if electrode == "positive": return "y" else: @@ -64,6 +135,11 @@ def _get_prefix(electrode): def _get_electrode_capacity_equation(options, electrode): + """ + Build equation for electrode capacity in composite electrodes. + + Returns Q = sum_i Q_i * (stoich_100_i - stoich_0_i) for all phases. + """ prefix = _get_prefix(electrode) e = electrode[0] i_am_composite = check_if_composite(options, electrode) @@ -86,6 +162,11 @@ def _get_electrode_capacity_equation(options, electrode): def _get_cyclable_lithium_equation(options, soc="100"): + """ + Build equation for total cyclable lithium in composite electrodes. + + Returns Q_Li = sum_i (Q_n_i * x_i + Q_p_i * y_i) for all phases at given SOC. + """ x_soc_1 = pybamm.Variable(f"x_{soc}_1", bounds=(0, 1)) y_soc_1 = pybamm.Variable(f"y_{soc}_1", bounds=(0, 1)) Q_n_1 = pybamm.InputParameter("Q_n_1") @@ -105,13 +186,109 @@ def _get_cyclable_lithium_equation(options, soc="100"): return lithium_primary_phases + lithium_secondary_phases +def _solve_secondary_stoichiometry( + primary_stoich: float, + parameter_values: pybamm.ParameterValues, + param: pybamm.LithiumIonParameters, + electrode: str, + direction: str | None, + options: dict, + T: float, + tol: float = 1e-6, +) -> float: + """ + Solve U_prim(z_1) = U_sec(z_2) to get z_2 given z_1. + + Parameters + ---------- + primary_stoich : float + The primary phase stoichiometry (x_1 or y_1) + parameter_values : pybamm.ParameterValues + The parameter values + param : pybamm.LithiumIonParameters + The parameter object + electrode : str + "negative" or "positive" + direction : str + "charge" or "discharge" + options : dict + Model options + T : float + Temperature + tol : float + Solver tolerance + + Returns + ------- + float + The secondary phase stoichiometry (x_2 or y_2) + """ + model = pybamm.BaseModel() + z_2 = pybamm.Variable("z_2", bounds=(0, 1)) + z_1 = pybamm.InputParameter("z_1") + + if electrode == "negative": + lith_prim = get_lithiation_delithiation( + direction, "negative", options, phase="primary" + ) + lith_sec = get_lithiation_delithiation( + direction, "negative", options, phase="secondary" + ) + U_prim = param.n.prim.U(z_1, T, lith_prim) + U_sec = param.n.sec.U(z_2, T, lith_sec) + else: + lith_prim = get_lithiation_delithiation( + direction, "positive", options, phase="primary" + ) + lith_sec = get_lithiation_delithiation( + direction, "positive", options, phase="secondary" + ) + U_prim = param.p.prim.U(z_1, T, lith_prim) + U_sec = param.p.sec.U(z_2, T, lith_sec) + + model.algebraic[z_2] = U_prim - U_sec + model.initial_conditions[z_2] = primary_stoich + model.variables["z_2"] = z_2 + + sim = pybamm.Simulation( + model, parameter_values=parameter_values, solver=pybamm.AlgebraicSolver(tol=tol) + ) + sol = sim.solve([0], inputs={"z_1": primary_stoich}) + return sol["z_2"].data[0] + + class ElectrodeSOHComposite(pybamm.BaseModel): """Model to calculate electrode-specific SOH for a cell with composite electrodes, adapted from :footcite:t:`Mohtat2019`. This model is mainly for internal use, to - calculate summary variables in a - simulation. + calculate summary variables in a simulation. Subscript 1 indicates primary phase and subscript 2 indicates secondary phase. + + The model calculates stoichiometries at three states: + + - 100% SOC (x_100, y_100): Limit state reached via charging + - 0% SOC (x_0, y_0): Limit state reached via discharging + - Initial SOC (x_init, y_init): Dynamic state, uses specified direction + + Stoichiometry limits (_100 and _0 variables) are evaluated in a + temperature-independent manner (at the reference temperature, ignoring entropy + effects), so that the SOC range is not temperature-dependent. In the presence of + a hysteresis model, equilibration in the charging direction is assumed for 100% + SOC (charging OCP branch for each material), and equilibration in the discharging + direction is assumed for 0% SOC (discharging OCP branch for each material). The + initial stoichiometries (_init variables) use the specified direction to account + for hysteresis during charge/discharge. + + Parameters + ---------- + options : dict + Model options including particle phases and OCP settings + direction : str, optional + "charge" or "discharge" - only affects initial stoichiometry calculation + name : str, optional + Model name (default: "ElectrodeSOH model") + initialization_method : str, optional + "voltage" or "SOC" (default: "voltage") """ def __init__( @@ -125,7 +302,6 @@ def __init__( super().__init__(name) param = pybamm.LithiumIonParameters(options) - # Start by just assuming known value is cyclable lithium (and solve all at once) Q_Li = pybamm.InputParameter("Q_Li") is_negative_composite = check_if_composite(options, "negative") is_positive_composite = check_if_composite(options, "positive") @@ -136,35 +312,45 @@ def __init__( y_0_1 = variables["y_0_1"] V_max = param.ocp_soc_100 V_min = param.ocp_soc_0 - # Here we use T_ref as the stoichiometry limits are defined using the reference - # state if is_negative_composite: x_100_2 = variables["x_100_2"] x_0_2 = variables["x_0_2"] self.algebraic[x_100_2] = param.n.sec.U( x_100_2, param.T_ref, - _get_lithiation_delithiation( - direction, "negative", options, phase="secondary" + get_lithiation_delithiation( + get_equilibrium_direction("100", "negative", options, "secondary"), + "negative", + options, + phase="secondary", ), ) - param.n.prim.U( x_100_1, param.T_ref, - _get_lithiation_delithiation( - direction, "negative", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("100", "negative", options, "primary"), + "negative", + options, + phase="primary", ), ) self.algebraic[x_0_2] = param.n.sec.U( x_0_2, param.T_ref, - _get_lithiation_delithiation( - direction, "negative", options, phase="secondary" + get_lithiation_delithiation( + get_equilibrium_direction("0", "negative", options, "secondary"), + "negative", + options, + phase="secondary", ), ) - param.n.prim.U( x_0_1, param.T_ref, - _get_lithiation_delithiation( - direction, "negative", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("0", "negative", options, "primary"), + "negative", + options, + phase="primary", ), ) if is_positive_composite: @@ -173,42 +359,60 @@ def __init__( self.algebraic[y_100_2] = param.p.sec.U( y_100_2, param.T_ref, - _get_lithiation_delithiation( - direction, "positive", options, phase="secondary" + get_lithiation_delithiation( + get_equilibrium_direction("100", "positive", options, "secondary"), + "positive", + options, + phase="secondary", ), ) - param.p.prim.U( y_100_1, param.T_ref, - _get_lithiation_delithiation( - direction, "positive", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("100", "positive", options, "primary"), + "positive", + options, + phase="primary", ), ) self.algebraic[y_0_2] = param.p.prim.U( y_0_1, param.T_ref, - _get_lithiation_delithiation( - direction, "positive", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("0", "positive", options, "primary"), + "positive", + options, + phase="primary", ), ) - param.p.sec.U( y_0_2, param.T_ref, - _get_lithiation_delithiation( - direction, "positive", options, phase="secondary" + get_lithiation_delithiation( + get_equilibrium_direction("0", "positive", options, "secondary"), + "positive", + options, + phase="secondary", ), ) self.algebraic[x_100_1] = ( param.p.prim.U( y_100_1, param.T_ref, - _get_lithiation_delithiation( - direction, "positive", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("100", "positive", options, "primary"), + "positive", + options, + phase="primary", ), ) - param.n.prim.U( x_100_1, param.T_ref, - _get_lithiation_delithiation( - direction, "negative", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("100", "negative", options, "primary"), + "negative", + options, + phase="primary", ), ) - V_max @@ -217,20 +421,25 @@ def __init__( param.p.prim.U( y_0_1, param.T_ref, - _get_lithiation_delithiation( - direction, "positive", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("0", "positive", options, "primary"), + "positive", + options, + phase="primary", ), ) - param.n.prim.U( x_0_1, param.T_ref, - _get_lithiation_delithiation( - direction, "negative", options, phase="primary" + get_lithiation_delithiation( + get_equilibrium_direction("0", "negative", options, "primary"), + "negative", + options, + phase="primary", ), ) - V_min ) - # arbitrary choice: use y_0_1 for the capacity equation self.algebraic[y_0_1] = _get_electrode_capacity_equation( options, "positive" ) - _get_electrode_capacity_equation(options, "negative") @@ -239,21 +448,19 @@ def __init__( x_init_1 = variables["x_init_1"] y_init_1 = variables["y_init_1"] if initialization_method == "voltage": - # Here we use T_init so that the initial voltage is correct, including the - # contribution from the entropic change V_init = pybamm.InputParameter("V_init") self.algebraic[x_init_1] = ( param.p.prim.U( y_init_1, param.T_init, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "positive", options, phase="primary" ), ) - param.n.prim.U( x_init_1, param.T_init, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "negative", options, phase="primary" ), ) @@ -284,23 +491,19 @@ def __init__( ) else: raise ValueError("Invalid initialization method") - # Add voltage equations for secondary phases (init), we use T_ref if setting - # based on SOC since the stoichiometry limits are defined using the reference - # state, and T_init if setting based on voltage since the entropic change is - # included in the voltage equation T = param.T_init if initialization_method == "voltage" else param.T_ref if is_positive_composite: y_init_2 = variables["y_init_2"] self.algebraic[y_init_2] = param.p.prim.U( y_init_1, T, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "positive", options, phase="primary" ), ) - param.p.sec.U( y_init_2, T, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "positive", options, phase="secondary" ), ) @@ -309,13 +512,13 @@ def __init__( self.algebraic[x_init_2] = param.n.prim.U( x_init_1, T, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "negative", options, phase="primary" ), ) - param.n.sec.U( x_init_2, T, - _get_lithiation_delithiation( + get_lithiation_delithiation( direction, "negative", options, phase="secondary" ), ) @@ -331,78 +534,444 @@ def __init__( def default_solver(self): return pybamm.AlgebraicSolver(method="lsq") + @staticmethod + def solve_split( + initial_value: float | str, + parameter_values: pybamm.ParameterValues, + direction: str | None = None, + param: pybamm.LithiumIonParameters | None = None, + options: dict | None = None, + tol: float = 1e-6, + inputs: dict | None = None, + ) -> dict[str, float]: + """ + Split solve approach for composite electrode SOH. -def get_initial_stoichiometries_composite( - initial_value, - parameter_values, - direction=None, - param=None, - options=None, - tol=1e-6, - inputs=None, - known_value="cyclable lithium capacity", - **kwargs, -): - """ - Get the minimum and maximum stoichiometries from the parameter values + Step 1: Solve for primary stoichiometries using non-composite model + Step 2: Solve U_prim(z_1) = U_sec(z_2) for secondary stoichiometries + + The equilibrium stoichiometries (x_100, x_0, y_100, y_0) are calculated + using direction=None (equilibrium branch). Only initial stoichiometries + (x_init, y_init) use the specified direction to account for hysteresis during + charge/discharge. + + Parameters + ---------- + initial_value : float or str + Target initial value. If float (0-1), interpreted as SOC. + If string ending in 'V', interpreted as voltage. + parameter_values : pybamm.ParameterValues + Parameter values for the simulation + direction : str, optional + "charge" or "discharge" for hysteresis direction (only affects + initial stoichiometries, not equilibrium values) + param : pybamm.LithiumIonParameters, optional + Parameter object + options : dict, optional + Model options + tol : float, optional + Solver tolerance (default 1e-6) + inputs : dict, optional + Additional inputs + + Returns + ------- + dict + Dictionary of stoichiometry values + """ + inputs = inputs or {} + param = param or pybamm.LithiumIonParameters(options) + + is_positive_composite = check_if_composite(options, "positive") + is_negative_composite = check_if_composite(options, "negative") + + Q_n_1 = parameter_values.evaluate(param.n.prim.Q_init, inputs=inputs) + Q_p_1 = parameter_values.evaluate(param.p.prim.Q_init, inputs=inputs) + Qs = {"Q_n_1": Q_n_1, "Q_p_1": Q_p_1} + if is_positive_composite: + Q_p_2 = parameter_values.evaluate(param.p.sec.Q_init, inputs=inputs) + Qs["Q_p_2"] = Q_p_2 + if is_negative_composite: + Q_n_2 = parameter_values.evaluate(param.n.sec.Q_init, inputs=inputs) + Qs["Q_n_2"] = Q_n_2 + + Q_Li = parameter_values.evaluate(param.Q_Li_particles_init, inputs=inputs) + + if isinstance(initial_value, str) and initial_value.endswith("V"): + V_init = float(initial_value[:-1]) + initialization_method = "voltage" + elif isinstance(initial_value, float) and 0 <= initial_value <= 1: + initialization_method = "SOC" + else: + raise ValueError( + "Invalid initial value. Expected a float between 0 and 1 " + "(for SOC) or a string ending in 'V' (for voltage), got " + f"{initial_value!r} of type {type(initial_value).__name__}" + ) + + Q_n_total = Q_n_1 + (Qs.get("Q_n_2", 0)) + Q_p_total = Q_p_1 + (Qs.get("Q_p_2", 0)) + + primary_options = _get_primary_only_options(options) + # _ElectrodeSOH uses get_equilibrium_direction internally for the equilibrium + # stoichiometries (x_0, x_100, y_0, y_100), consistent with the full solve + primary_model = _ElectrodeSOH( + direction=direction, + param=param, + solve_for=["x_0", "x_100"], + known_value="cyclable lithium capacity", + options=primary_options, + ) + + primary_inputs = {**inputs, "Q_n": Q_n_total, "Q_p": Q_p_total, "Q_Li": Q_Li} + primary_sim = pybamm.Simulation( + primary_model, + parameter_values=parameter_values, + solver=pybamm.AlgebraicSolver(tol=tol), + ) + + primary_sol = primary_sim.solve([0], inputs=primary_inputs) + x_100_1 = primary_sol["x_100"].data[0] + x_0_1 = primary_sol["x_0"].data[0] + y_100_1 = primary_sol["y_100"].data[0] + y_0_1 = primary_sol["y_0"].data[0] + + T_ref = parameter_values["Reference temperature [K]"] + result = { + "x_100_1": x_100_1, + "x_0_1": x_0_1, + "y_100_1": y_100_1, + "y_0_1": y_0_1, + } + + if is_negative_composite: + result["x_100_2"] = _solve_secondary_stoichiometry( + x_100_1, + parameter_values, + param, + "negative", + "charge", # 100% SOC is reached via charging + options, + T_ref, + tol, + ) + result["x_0_2"] = _solve_secondary_stoichiometry( + x_0_1, + parameter_values, + param, + "negative", + "discharge", # 0% SOC is reached via discharging + options, + T_ref, + tol, + ) + + if is_positive_composite: + result["y_100_2"] = _solve_secondary_stoichiometry( + y_100_1, + parameter_values, + param, + "positive", + "charge", # 100% SOC is reached via charging + options, + T_ref, + tol, + ) + result["y_0_2"] = _solve_secondary_stoichiometry( + y_0_1, + parameter_values, + param, + "positive", + "discharge", # 0% SOC is reached via discharging + options, + T_ref, + tol, + ) + + if initialization_method == "voltage": + T_init = parameter_values["Initial temperature [K]"] + soc_model = pybamm.BaseModel() + x_init = pybamm.Variable("x_init", bounds=(0, 1)) + y_init = y_0_1 + (x_init - x_0_1) / (x_100_1 - x_0_1) * (y_100_1 - y_0_1) + lith_pos = get_lithiation_delithiation( + direction, "positive", options, phase="primary" + ) + lith_neg = get_lithiation_delithiation( + direction, "negative", options, phase="primary" + ) + Up = param.p.prim.U(y_init, T_init, lith_pos) + Un = param.n.prim.U(x_init, T_init, lith_neg) + soc_model.algebraic[x_init] = Up - Un - V_init + soc_model.initial_conditions[x_init] = (x_0_1 + x_100_1) / 2 + soc_model.variables["x_init"] = x_init + soc_model.variables["y_init"] = y_init + soc_sim = pybamm.Simulation( + soc_model, + parameter_values=parameter_values, + solver=pybamm.AlgebraicSolver(tol=tol), + ) + soc_sol = soc_sim.solve([0], inputs=inputs) + x_init_1 = soc_sol["x_init"].data[0] + y_init_1 = soc_sol["y_init"].data[0] + else: # SOC initialization + soc = initial_value + x_init_1 = x_0_1 + soc * (x_100_1 - x_0_1) + y_init_1 = y_0_1 + soc * (y_100_1 - y_0_1) + + result["x_init_1"] = x_init_1 + result["y_init_1"] = y_init_1 + + T_for_init = ( + parameter_values["Initial temperature [K]"] + if initialization_method == "voltage" + else T_ref + ) + if is_negative_composite: + result["x_init_2"] = _solve_secondary_stoichiometry( + x_init_1, + parameter_values, + param, + "negative", + direction, + options, + T_for_init, + tol, + ) + if is_positive_composite: + result["y_init_2"] = _solve_secondary_stoichiometry( + y_init_1, + parameter_values, + param, + "positive", + direction, + options, + T_for_init, + tol, + ) + + return result + + @staticmethod + def solve_full( + initial_value: float | str, + parameter_values: pybamm.ParameterValues, + direction: str | None = None, + param: pybamm.LithiumIonParameters | None = None, + options: dict | None = None, + tol: float = 1e-6, + inputs: dict | None = None, + initial_conditions: dict[str, float] | None = None, + ) -> dict[str, float]: + """ + Full solve approach: solve all stoichiometries simultaneously. + + Uses the full ElectrodeSOHComposite algebraic model to solve for all + stoichiometries at once. The equilibrium stoichiometries (x_100, x_0, + y_100, y_0) are calculated using direction=None (equilibrium branch). + Only initial stoichiometries (x_init, y_init) use the specified direction to + account for hysteresis. Parameters ---------- - initial_value : float - Target initial value. - If integer, interpreted as SOC, must be between 0 and 1. - If string e.g. "4 V", interpreted as voltage, - must be between V_min and V_max. + initial_value : float or str + Target initial value. If float (0-1), interpreted as SOC. + If string ending in 'V', interpreted as voltage. + parameter_values : pybamm.ParameterValues + Parameter values for the simulation direction : str, optional - The OCV branch to use in the electrode SOH model. Can be "charge" or - "discharge". + "charge" or "discharge" for hysteresis direction (only affects + initial stoichiometries, not equilibrium values) + param : pybamm.LithiumIonParameters, optional + Parameter object + options : dict, optional + Model options tol : float, optional - The tolerance for the solver used to compute the initial stoichiometries. - A lower value results in higher precision but may increase computation time. - Default is 1e-6. + Solver tolerance (default 1e-6) inputs : dict, optional - A dictionary of input parameters passed to the model. - known_value : str, optional - The known value needed to complete the electrode SOH model. - Can be "cyclable lithium capacity". + Additional inputs + initial_conditions : dict, optional + Dictionary of initial conditions for variables (e.g., from split solve) + + Returns + ------- + dict + Dictionary of stoichiometry values + """ + inputs = inputs or {} + param = param or pybamm.LithiumIonParameters(options) + + is_positive_composite = check_if_composite(options, "positive") + is_negative_composite = check_if_composite(options, "negative") + + Q_n_1 = parameter_values.evaluate(param.n.prim.Q_init, inputs=inputs) + Q_p_1 = parameter_values.evaluate(param.p.prim.Q_init, inputs=inputs) + Qs = {"Q_n_1": Q_n_1, "Q_p_1": Q_p_1} + if is_positive_composite: + Q_p_2 = parameter_values.evaluate(param.p.sec.Q_init, inputs=inputs) + Qs["Q_p_2"] = Q_p_2 + if is_negative_composite: + Q_n_2 = parameter_values.evaluate(param.n.sec.Q_init, inputs=inputs) + Qs["Q_n_2"] = Q_n_2 + + Q_Li = parameter_values.evaluate(param.Q_Li_particles_init, inputs=inputs) + + if isinstance(initial_value, str) and initial_value.endswith("V"): + V_init = float(initial_value[:-1]) + initialization_method = "voltage" + elif isinstance(initial_value, float) and 0 <= initial_value <= 1: + initialization_method = "SOC" + else: + raise ValueError( + "Invalid initial value. Expected a float between 0 and 1 " + "(for SOC) or a string ending in 'V' (for voltage), got " + f"{initial_value!r} of type {type(initial_value).__name__}" + ) + + all_inputs = {**inputs, **Qs, "Q_Li": Q_Li} + if initialization_method == "voltage": + all_inputs["V_init"] = V_init + else: + all_inputs["SOC_init"] = initial_value + + model = ElectrodeSOHComposite( + options, direction, initialization_method=initialization_method + ) + sim = pybamm.Simulation( + model, + parameter_values=parameter_values, + solver=pybamm.AlgebraicSolver(tol=tol), + ) + + if initial_conditions is not None: + sim.build() + sim.built_model.set_initial_conditions_from( + initial_conditions, inputs=all_inputs + ) + sol = sim.solve([0], inputs=all_inputs) + + return {var: sol[var].entries[0] for var in model.variables.keys()} + + +def get_initial_stoichiometries_composite( + initial_value: float | str, + parameter_values: pybamm.ParameterValues, + direction: str | None = None, + param: pybamm.LithiumIonParameters | None = None, + options: dict | None = None, + tol: float = 1e-6, + inputs: dict | None = None, + known_value: str = "cyclable lithium capacity", + try_split_solve: bool = True, + **kwargs: Any, +) -> dict[str, float]: + """ + Get the stoichiometries for composite electrodes from parameter values. + + Calculates stoichiometries at three states: + - 100% SOC (x_100, y_100): Equilibrium state (direction=None) + - 0% SOC (x_0, y_0): Equilibrium state (direction=None) + - Initial SOC (x_init, y_init): Dynamic state (uses specified direction) + + For electrode models with OCP hysteresis, the equilibrium stoichiometries use + the charging OCP branch at 100% SOC and the discharging OCP branch at 0% SOC. + For models without hysteresis there is only one OCP curve, so the direction + does not affect the equilibrium stoichiometries. + + Parameters + ---------- + initial_value : float or str + Target initial value. + If float between 0 and 1, interpreted as SOC. + If string ending in 'V' (e.g., "4 V"), interpreted as voltage, + must be between V_min and V_max. + parameter_values : pybamm.ParameterValues + Parameter values for the simulation + direction : str, optional + The OCV branch to use for initial stoichiometries. Can be "charge" or + "discharge". Only affects x_init/y_init, not equilibrium values. + param : pybamm.LithiumIonParameters, optional + Parameter object + options : dict, optional + Model options + tol : float, optional + The tolerance for the solver used to compute the initial stoichiometries. + A lower value results in higher precision but may increase computation time. + Default is 1e-6. + inputs : dict, optional + A dictionary of input parameters passed to the model. + known_value : str, optional + The known value needed to complete the electrode SOH model. + Can be "cyclable lithium capacity". + try_split_solve : bool, optional + Whether to use the split solve method to improve robustness. Default is True. + When True, if the full solve fails: + 1. Run split solve to get approximate stoichiometries + 2. Use these as initial conditions and retry full solve + 3. If retry succeeds, return full solve results + 4. If retry fails, raise error + + Returns + ------- + dict + Dictionary of stoichiometry values for all phases at 0%, 100%, and initial SOC """ inputs = inputs or {} + param = param or pybamm.LithiumIonParameters(options) + if known_value != "cyclable lithium capacity": raise ValueError( "Only `cyclable lithium capacity` is supported for composite electrodes" ) - Q_n_1 = parameter_values.evaluate(param.n.prim.Q_init, inputs=inputs) - Q_p_1 = parameter_values.evaluate(param.p.prim.Q_init, inputs=inputs) - is_positive_composite = check_if_composite(options, "positive") - is_negative_composite = check_if_composite(options, "negative") - Qs = { - "Q_n_1": Q_n_1, - "Q_p_1": Q_p_1, - } - if is_positive_composite: - Q_p_2 = parameter_values.evaluate(param.p.sec.Q_init, inputs=inputs) - Qs["Q_p_2"] = Q_p_2 - if is_negative_composite: - Q_n_2 = parameter_values.evaluate(param.n.sec.Q_init, inputs=inputs) - Qs["Q_n_2"] = Q_n_2 - - Q_Li = parameter_values.evaluate(param.Q_Li_particles_init, inputs=inputs) - all_inputs = {**inputs, **Qs, "Q_Li": Q_Li} - # Solve the model and check outputs - if isinstance(initial_value, str) and initial_value.endswith("V"): - all_inputs["V_init"] = float(initial_value[:-1]) - initialization_method = "voltage" - elif isinstance(initial_value, float) and initial_value >= 0 and initial_value <= 1: - initialization_method = "SOC" - all_inputs["SOC_init"] = initial_value - else: - raise ValueError("Invalid initial value") - model = ElectrodeSOHComposite( - options, direction, initialization_method=initialization_method - ) - sim = pybamm.Simulation( - model, parameter_values=parameter_values, solver=pybamm.AlgebraicSolver(tol=tol) - ) - sol = sim.solve([0, 1], inputs=all_inputs) - return {var: sol[var].entries[0] for var in model.variables.keys()} + try: + return ElectrodeSOHComposite.solve_full( + initial_value, + parameter_values, + direction=direction, + param=param, + options=options, + tol=tol, + inputs=inputs, + ) + except (pybamm.SolverError, ValueError) as first_error: + if try_split_solve: + try: + split_results = ElectrodeSOHComposite.solve_split( + initial_value, + parameter_values, + direction=direction, + param=param, + options=options, + tol=tol, + inputs=inputs, + ) + + try: + return ElectrodeSOHComposite.solve_full( + initial_value, + parameter_values, + direction=direction, + param=param, + options=options, + tol=tol, + inputs=inputs, + initial_conditions=split_results, + ) + except (pybamm.SolverError, ValueError) as retry_error: + raise ValueError( + f"Failed to solve composite electrode SOH. " + f"Initial full solve error: {first_error}. " + f"Retry with split solve initial conditions also failed: " + f"{retry_error}" + ) from retry_error + + except (pybamm.SolverError, ValueError) as split_error: + raise ValueError( + f"Failed to solve composite electrode SOH. " + f"Full solve error: {first_error}. " + f"Split solve error: {split_error}" + ) from split_error + else: + raise ValueError( + f"Failed to solve composite electrode SOH: {first_error}" + ) from first_error diff --git a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py index 10e3185f40..19c6330313 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py @@ -3,7 +3,7 @@ # import pybamm -from .util import _get_lithiation_delithiation, check_if_composite +from .util import check_if_composite, get_lithiation_delithiation class ElectrodeSOHHalfCell(pybamm.BaseModel): @@ -37,7 +37,7 @@ def __init__(self, name="ElectrodeSOH model", direction=None, options=None): Q_w = pybamm.InputParameter("Q_w") T_ref = param.T_ref U_w = param.p.prim.U - lith_delith_primary = _get_lithiation_delithiation( + lith_delith_primary = get_lithiation_delithiation( direction, "positive", options, phase="primary" ) @@ -49,7 +49,7 @@ def __init__(self, name="ElectrodeSOH model", direction=None, options=None): Q_w_2 = pybamm.InputParameter("Q_w_2") U_w_2 = param.p.sec.U Q_2 = Q_w_2 * (x_0_2 - x_100_2) - lith_delith_secondary = _get_lithiation_delithiation( + lith_delith_secondary = get_lithiation_delithiation( direction, "positive", options, phase="secondary" ) @@ -141,12 +141,12 @@ def get_initial_stoichiometry_half_cell( ) x_0, x_100 = x_dict["x_0"], x_dict["x_100"] is_composite = check_if_composite(options, "positive") - lith_delith_primary = _get_lithiation_delithiation( + lith_delith_primary = get_lithiation_delithiation( direction, "positive", options, phase="primary" ) if is_composite: x_0_2, x_100_2 = x_dict["x_0_2"], x_dict["x_100_2"] - lith_delith_secondary = _get_lithiation_delithiation( + lith_delith_secondary = get_lithiation_delithiation( direction, "positive", options, phase="secondary" ) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/util.py b/src/pybamm/models/full_battery_models/lithium_ion/util.py index 545a8e055a..c614e2ea15 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/util.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/util.py @@ -38,12 +38,30 @@ def _has_hysteresis(electrode, options, phase=None): for item in ocp_option ) else: - raise ValueError( - f"Invalid open-circuit potential option: {domain_options['open-circuit potential']}" - ) + ocp_opt = domain_options["open-circuit potential"] + raise ValueError(f"Invalid open-circuit potential option: {ocp_opt}") -def _get_lithiation_delithiation(direction, electrode, options, phase=None): +def get_lithiation_delithiation(direction, electrode, options, phase=None): + """ + Get the lithiation/delithiation direction for OCP evaluation. + + Parameters + ---------- + direction : str or None + The cell-level direction: "charge", "discharge", or None (equilibrium). + electrode : str + "negative" or "positive" + options : dict or BatteryModelOptions + Model options containing OCP settings + phase : str, optional + "primary" or "secondary" for composite electrodes + + Returns + ------- + str or None + "lithiation", "delithiation", or None (for equilibrium/no hysteresis) + """ if direction is None or not _has_hysteresis(electrode, options, phase): return None elif (direction == "charge" and electrode == "negative") or ( @@ -56,3 +74,44 @@ def _get_lithiation_delithiation(direction, electrode, options, phase=None): return "delithiation" else: raise ValueError() + + +def get_equilibrium_direction(soc_state, electrode, options, phase=None): + """ + Get the appropriate cell direction for equilibrium stoichiometry calculation. + + For electrodes with hysteresis, equilibrium stoichiometries should be calculated + on the OCP branch corresponding to how the cell reached that SOC: + - 100% SOC: reached via charging → use "charge" direction + - 0% SOC: reached via discharging → use "discharge" direction + + For electrodes without hysteresis, returns None (use equilibrium OCP). + + Parameters + ---------- + soc_state : str + "100" for 100% SOC or "0" for 0% SOC + electrode : str + "negative" or "positive" + options : dict or BatteryModelOptions + Model options containing OCP settings + phase : str, optional + "primary" or "secondary" for composite electrodes + + Returns + ------- + str or None + "charge" (for 100% SOC with hysteresis), "discharge" (for 0% SOC with + hysteresis), or None (no hysteresis) + """ + if not _has_hysteresis(electrode, options, phase): + return None + + if soc_state == "100": + # 100% SOC is reached by charging + return "charge" + elif soc_state == "0": + # 0% SOC is reached by discharging + return "discharge" + else: + raise ValueError(f"soc_state must be '100' or '0', got '{soc_state}'") diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py index 894091e997..1229b7d130 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py @@ -248,7 +248,7 @@ def _get_params_and_options(composite_electrode): options = {"particle phases": phases} return params, options - @pytest.mark.parametrize("initial_value", ["4.0 V", 0.5]) + @pytest.mark.parametrize("initial_value", ["2.7V", "4.0 V", 0.1, 0.5, 0.9]) @pytest.mark.parametrize( "composite_electrode", [ @@ -285,7 +285,8 @@ def test_half_cell_with_same_ocp_curves(self, composite_electrode, initial_value - param.n.prim.U(results["x_init_1"], param.T_ref) ) == pytest.approx(4.0, abs=1e-05) - def test_chen2020_composite_defaults(self): + @pytest.mark.parametrize("initial_value", ["2.7V", "4.0 V", 0.1, 0.5, 0.9]) + def test_chen2020_composite_defaults(self, initial_value): pvals = pybamm.ParameterValues("Chen2020_composite") options = {"particle phases": ("2", "1")} param = pybamm.LithiumIonParameters(options=options) @@ -298,22 +299,77 @@ def test_chen2020_composite_defaults(self): } ) results = pybamm.lithium_ion.get_initial_stoichiometries_composite( - "4.0 V", pvals, param=param, options=options, tol=1e-6, direction=None + initial_value, pvals, param=param, options=options, tol=1e-6, direction=None ) # Basic sanity: solution includes expected variables and bounded stoichiometries for key, val in results.items(): if key.startswith(("x_", "y_")): assert 0 <= val <= 1 pvals_set = pybamm.lithium_ion.set_initial_state( - "4.0 V", pvals, param=param, options=options, tol=1e-6 + initial_value, pvals, param=param, options=options, tol=1e-6 ) - assert pvals_set.evaluate( - param.p.prim.U(results["y_init_1"], param.T_ref) - - param.n.prim.U(results["x_init_1"], param.T_ref) - ) == pytest.approx(4.0, abs=1e-05) + if isinstance(initial_value, str): + V_target = float(initial_value[:-1]) + assert pvals_set.evaluate( + param.p.prim.U(results["y_init_1"], param.T_ref) + - param.n.prim.U(results["x_init_1"], param.T_ref) + ) == pytest.approx(V_target, abs=1e-05) + + def test_chen2020_composite_default_solve(self): + pvals = pybamm.ParameterValues("Chen2020_composite") + options = {"particle phases": ("2", "1")} + param = pybamm.LithiumIonParameters(options=options) + pvals.update( + { + "Secondary: Initial concentration in negative electrode [mol.m-3]": 2.3512e05 + } + ) + + # Test default solve with split fallback enabled + results_with_fallback = ( + pybamm.lithium_ion.get_initial_stoichiometries_composite( + "4.0 V", + pvals, + param=param, + options=options, + tol=1e-6, + direction=None, + try_split_solve=True, + ) + ) + + # Test without fallback (should also work for this case) + results_no_fallback = pybamm.lithium_ion.get_initial_stoichiometries_composite( + "4.0 V", + pvals, + param=param, + options=options, + tol=1e-6, + direction=None, + try_split_solve=False, + ) + + # Both should give similar results for this well-behaved case + for key in results_with_fallback: + if key.startswith(("x_", "y_")): + assert results_with_fallback[key] == pytest.approx( + results_no_fallback[key], abs=1e-5 + ), ( + f"Mismatch for {key}: " + f"with_fallback={results_with_fallback[key]}, " + f"no_fallback={results_no_fallback[key]}" + ) def test_chen2020_composite_defaults_hysteresis(self): pvals = pybamm.ParameterValues("Chen2020_composite") + # Solving ESOH with the original Chen2020_composite parameters gives a 0% SOC + # voltage of 2.53V not 2.5V. We fix this by reducing the secondary initial + # concentration, which adjusts Q_Li to make the system consistent. + pvals.update( + { + "Secondary: Initial concentration in negative electrode [mol.m-3]": 2.3512e05 + } + ) options = { "particle phases": ("2", "1"), "open-circuit potential": (("single", "current sigmoid"), "single"), @@ -324,17 +380,55 @@ def test_chen2020_composite_defaults_hysteresis(self): pvals, param=param, options=options, - tol=1e-1, + tol=1e-6, direction="discharge", ) results_charge = pybamm.lithium_ion.get_initial_stoichiometries_composite( - "4.0 V", pvals, param=param, options=options, tol=1e-1, direction="charge" + "4.0 V", pvals, param=param, options=options, tol=1e-6, direction="charge" ) # Basic sanity: solution includes expected variables and bounded stoichiometries - for key, val in results_discharge.items(): + for key in results_discharge: if key.startswith(("x_", "y_")): - assert 0 <= val <= 1 - assert results_discharge[key] != results_charge[key] + val_discharge = results_discharge[key] + val_charge = results_charge[key] + # Check both values are bounded + assert 0 <= val_discharge <= 1 + assert 0 <= val_charge <= 1 + # Check values are approximately equal if not 'init' in key (same + # voltage may be different stoichiometries due to hysteresis) + if "init" not in key: + assert val_discharge == pytest.approx(val_charge, abs=1e-04) + + def test_chen2020_composite_solve_with_hysteresis(self): + pvals = pybamm.ParameterValues("Chen2020_composite") + pvals.update( + { + "Secondary: Initial concentration in negative electrode [mol.m-3]": 2.3512e05 + } + ) + options = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + } + param = pybamm.LithiumIonParameters(options=options) + + # Test that solver works with hysteresis in both directions + for direction in ["charge", "discharge"]: + results = pybamm.lithium_ion.get_initial_stoichiometries_composite( + "4.0 V", + pvals, + param=param, + options=options, + tol=1e-6, + direction=direction, + ) + + # Verify all stoichiometries are within bounds + for key in results: + if key.startswith(("x_", "y_")): + assert 0 <= results[key] <= 1, ( + f"Stoichiometry {key} out of bounds: {results[key]}" + ) class TestElectrodeSOHMSMR: From 2641f1c728afdb2547ace0f484ba37545e0196d5 Mon Sep 17 00:00:00 2001 From: bradyplanden Date: Wed, 21 Jan 2026 09:42:09 +0000 Subject: [PATCH 3/3] Increment version to v25.12.1 --- CHANGELOG.md | 7 +++++++ CITATION.cff | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2637bbb731..f1b0f573eb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) +## Features +- Improve the performance of matrix multiplication with CasADi expressions. ([#5351](https://github.com/pybamm-team/PyBaMM/pull/5351)) + +## Bug fixes + +# [v25.12.1](https://github.com/pybamm-team/PyBaMM/tree/v25.12.1) - 2026-01-21 + ## Features - Improve the robustness of the `ElectrodeSOHComposite` model and ensure the calculation of equilibrium stoichiometries is consistent for models with hysteresis. ([#5347](https://github.com/pybamm-team/PyBaMM/pull/5347)) diff --git a/CITATION.cff b/CITATION.cff index f58eefad5c..56445f4567 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -24,6 +24,6 @@ keywords: - "expression tree" - "python" - "symbolic differentiation" -version: "25.12.0" +version: "25.12.1" repository-code: "https://github.com/pybamm-team/PyBaMM" title: "Python Battery Mathematical Modelling (PyBaMM)"