Repository: pybamm-team/PyBaMM
Pull request #4983 (2025-04-17 20:24:17)

Release v25.4.2
Bug fix release

Note: Must be merged not squashed


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0ffdcc782..15b079848 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,10 +1,20 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+# [v25.4.2](https://github.com/pybamm-team/PyBaMM/tree/v25.4.2) - 2025-04-17
+
+## Bug fixes
+
+- Improve reliability of `AlgebraicSolver` and change `ElectrodeSOHHalfCell` solver to a Trust-Region method. ([#4982](https://github.com/pybamm-team/PyBaMM/pull/4982))
+
 # [v25.4.1](https://github.com/pybamm-team/PyBaMM/tree/v25.4.1) - 2025-04-16
 
+## Features
+
+- Added "use lumped thermal capacity" option in lumped thermal model ([#4968](https://github.com/pybamm-team/PyBaMM/pull/4968))
+
 ## Bug fixes
 
-- Remove a regularization term in the harmonic mean.  ([#4977](https://github.com/pybamm-team/PyBaMM/pull/4977))
+- Remove a regularization term in the harmonic mean. ([#4977](https://github.com/pybamm-team/PyBaMM/pull/4977))
 
 ## Breaking changes
 
@@ -64,7 +74,6 @@
 - Added `Symbol.post_order()` method to return an iterable that steps through the tree in post-order fashion. ([#4684](https://github.com/pybamm-team/PyBaMM/pull/4684))
 - Porosity change now works for composite electrode ([#4417](https://github.com/pybamm-team/PyBaMM/pull/4417))
 - Added two more submodels (options) for the SEI: Lars von Kolzenberg (2020) model and Tunneling Limit model ([#4394](https://github.com/pybamm-team/PyBaMM/pull/4394))
-- Added "use lumped thermal capacity" option in lumped thermal model ([#4968](https://github.com/pybamm-team/PyBaMM/pull/4968))
 
 ## Breaking changes
 
diff --git a/CITATION.cff b/CITATION.cff
index 6fdc6e550..4c6784413 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -24,6 +24,6 @@ keywords:
   - "expression tree"
   - "python"
   - "symbolic differentiation"
-version: "25.4.1"
+version: "25.4.2"
 repository-code: "https://github.com/pybamm-team/PyBaMM"
 title: "Python Battery Mathematical Modelling (PyBaMM)"
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 012169ef0..ef29feeb7 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
@@ -54,7 +54,7 @@ class ElectrodeSOHHalfCell(pybamm.BaseModel):
     @property
     def default_solver(self):
         # Use AlgebraicSolver as CasadiAlgebraicSolver gives unnecessary warnings
-        return pybamm.AlgebraicSolver(method="minimize  L-BFGS-B", tol=1e-7)
+        return pybamm.AlgebraicSolver(method="lsq__trf", tol=1e-7)
 
 
 def get_initial_stoichiometry_half_cell(
diff --git a/src/pybamm/solvers/algebraic_solver.py b/src/pybamm/solvers/algebraic_solver.py
index 9b6663d00..8a6755c18 100644
--- a/src/pybamm/solvers/algebraic_solver.py
+++ b/src/pybamm/solvers/algebraic_solver.py
@@ -28,13 +28,29 @@ class AlgebraicSolver(pybamm.BaseSolver):
         Any options to pass to the rootfinder. Vary depending on which method is chosen.
         Please consult `SciPy documentation
         <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.show_options.html>`_
-        for details.
+        for details. Addititional options to pass to the solver, by default:
+
+        .. code-block:: python
+
+            extra_options = {
+                # Tolerance for termination by the change of the independent variables.
+                "xtol": 1e-12,
+                # Tolerance for termination by the norm of the gradient.
+                "gtol": 1e-12,
+            }
     """
 
     def __init__(self, method="lm", tol=1e-6, extra_options=None):
         super().__init__(method=method)
         self.tol = tol
-        self.extra_options = extra_options or {}
+
+        default_extra_options = {
+            "xtol": 1e-12,
+            "gtol": 1e-12,
+        }
+        extra_options = extra_options or {}
+        self.extra_options = default_extra_options | extra_options
+
         self.name = f"Algebraic solver ({method})"
         self._algebraic_solver = True
         pybamm.citations.register("Virtanen2020")
@@ -158,6 +174,7 @@ class AlgebraicSolver(pybamm.BaseSolver):
                         **self.extra_options,
                     )
                     integration_time += timer.time()
+                    success |= sol.success
                 # Methods which use minimize are specified as either "minimize",
                 # which uses the default method, or with "minimize__methodname"
                 elif self.method.startswith("minimize"):
@@ -183,7 +200,14 @@ class AlgebraicSolver(pybamm.BaseSolver):
                         bounds = [
                             (lb, ub) for lb, ub in zip(model.bounds[0], model.bounds[1])
                         ]
-                        extra_options["bounds"] = bounds
+                    else:
+                        bounds = None
+
+                    hess = extra_options.get("hess", None)
+                    hessp = extra_options.get("hessp", None)
+                    constraints = extra_options.get("constraints", ())
+                    callback = extra_options.get("callback", None)
+
                     timer.reset()
                     sol = optimize.minimize(
                         root_norm,
@@ -191,9 +215,16 @@ class AlgebraicSolver(pybamm.BaseSolver):
                         method=method,
                         tol=self.tol,
                         jac=jac_norm,
-                        **extra_options,
+                        bounds=bounds,
+                        hess=hess,
+                        hessp=hessp,
+                        constraints=constraints,
+                        callback=callback,
+                        options=extra_options,
                     )
                     integration_time += timer.time()
+                    # The solution.success flag is unreliable, so we ignore
+                    # it and use the norm of the residual instead
                 else:
                     timer.reset()
                     sol = optimize.root(
@@ -205,25 +236,21 @@ class AlgebraicSolver(pybamm.BaseSolver):
                         options=self.extra_options,
                     )
                     integration_time += timer.time()
+                    # The solution.success flag is unreliable, so we ignore
+                    # it and use the norm of the residual instead
 
-                if sol.success and np.all(abs(sol.fun) < self.tol):
-                    # update initial guess for the next iteration
-                    y0_alg = sol.x
-                    # update solution array
+                y0_alg = sol.x
+                success |= np.all(abs(sol.fun) < self.tol)
+                if success:
                     y_alg[:, idx] = y0_alg
-                    success = True
-                elif not sol.success:
+                    break
+
+                if itr > maxiter:
                     raise pybamm.SolverError(
-                        f"Could not find acceptable solution: {sol.message}"
+                        "Could not find acceptable solution: solver terminated "
+                        "unsuccessfully and maximum solution error "
+                        f"({np.max(abs(sol.fun))}) above tolerance ({self.tol})"
                     )
-                else:
-                    y0_alg = sol.x
-                    if itr > maxiter:
-                        raise pybamm.SolverError(
-                            "Could not find acceptable solution: solver terminated "
-                            "successfully, but maximum solution error "
-                            f"({np.max(abs(sol.fun))}) above tolerance ({self.tol})"
-                        )
                 itr += 1
 
         # Concatenate differential part


Test patch:
diff --git a/tests/unit/test_solvers/test_algebraic_solver.py b/tests/unit/test_solvers/test_algebraic_solver.py
index 556b4de90..685e39833 100644
--- a/tests/unit/test_solvers/test_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_algebraic_solver.py
@@ -14,7 +14,7 @@ class TestAlgebraicSolver:
             method="hybr", tol=1e-4, extra_options={"maxfev": 100}
         )
         assert solver.method == "hybr"
-        assert solver.extra_options == {"maxfev": 100}
+        assert solver.extra_options == {"xtol": 1e-12, "gtol": 1e-12, "maxfev": 100}
         assert solver.tol == 1e-4
 
         solver.method = "krylov"
@@ -78,14 +78,14 @@ class TestAlgebraicSolver:
         solver = pybamm.AlgebraicSolver(method="hybr")
         with pytest.raises(
             pybamm.SolverError,
-            match="Could not find acceptable solution: The iteration is not making",
+            match="Could not find acceptable solution",
         ):
             solver._integrate(model, np.array([0]))
 
         solver = pybamm.AlgebraicSolver()
         with pytest.raises(
             pybamm.SolverError,
-            match="Could not find acceptable solution: solver terminated",
+            match="Could not find acceptable solution",
         ):
             solver._integrate(model, np.array([0]))
 
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index f67248ade..ae407e422 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -247,13 +247,13 @@ class TestBaseSolver:
 
         with pytest.raises(
             pybamm.SolverError,
-            match="Could not find acceptable solution: The iteration is not making",
+            match="Could not find acceptable solution",
         ):
             solver.calculate_consistent_state(Model())
         solver = pybamm.BaseSolver(root_method="lm")
         with pytest.raises(
             pybamm.SolverError,
-            match="Could not find acceptable solution: solver terminated",
+            match="Could not find acceptable solution",
         ):
             solver.calculate_consistent_state(Model())
         # with casadi