Repository: pybamm-team/PyBaMM
Pull request #4982 (2025-04-17 19:27:21)

Update `AlgebraicSolver` tolerances + default `ESOH` solver method
# Description

Previously, the `AlgebraicSolver` would only declare success if the absolute tolerance of the residuals were below the `tol` threshold. However, the scipy solvers would frequently terminate if other thresholds were hit, such as `xtol` or `gtol`, and `AlgebraicSolver` may erroneously throw an error even after a successful solve

This PR:
- tightens `xtol` and `gtol` from `1e-8` to `1e-12` to reflect the fact that we desire `all(abs(f(x)) < tol)`
- updates the success criterion to terminate either if the solver successfully exits or `all(abs(f(x)) < tol)` is satisfied
- also changes the default solver in `ElectrodeSOHHalfCell` to trust-region instead of a minimizer

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #)

# Important checks:

Please confirm the following before marking the PR as ready for review:
- No style issues: `nox -s pre-commit`
- All tests pass: `nox -s tests`
- The documentation builds: `nox -s doctests`
- Code is commented for hard-to-understand areas
- Tests added that prove fix is effective or that feature works


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4982?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 98.62%. Comparing base [(`4363962`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/4363962cba8ccf5028227a27922393ba8b9a827a?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`5f028da`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/5f028da2bdc1e4365258f5e588786994c0b957ec?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #4982      +/-   ##
===========================================
- Coverage    98.64%   98.62%   -0.03%     
===========================================
  Files          304      304              
  Lines        23750    23755       +5     
===========================================
- Hits         23429    23428       -1     
- Misses         321      327       +6     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4982?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

<details><summary> :rocket: New features to boost your workflow: </summary>

- :snowflake: [Test Analytics](https://docs.codecov.com/docs/test-analytics): Detect flaky tests, report on failures, and find test suite problems.
</details>

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0ffdcc782..5f2979061 100644
--- a/CHANGELOG.md
+++ b/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
 
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