Update AlgebraicSolver tolerances + default ESOH solver method (#4982)

* update success criterion + tighter tols

* tests

* change esoh from minimizer to root finder

* docs

* Update CHANGELOG.md

Author: MarcBerliner

CHANGELOG.md

@@ -1,10 +1,14 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Bug fixes
+
+- Improve reliablity 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
 
 ## 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
 

src/pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py

@@ -54,7 +54,7 @@ def __init__(self, name="ElectrodeSOH model"):
     @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(

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 @@ def jac_fn(y_alg, jac=jac):
                         **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,17 +200,31 @@ def jac_norm(y, jac_fn=jac_fn):
                         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,
                         y0_alg,
                         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 @@ def jac_norm(y, jac_fn=jac_fn):
                         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

tests/unit/test_solvers/test_algebraic_solver.py

@@ -14,7 +14,7 @@ def test_algebraic_solver_init(self):
             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 @@ def algebraic_eval(self, t, y, inputs):
         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]))
 

tests/unit/test_solvers/test_base_solver.py

@@ -247,13 +247,13 @@ def algebraic_eval(self, t, y, inputs):
 
         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