From 53a8040c79cc6990169a34ea2779f260e049b10e Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 11 Oct 2022 18:35:31 -0400 Subject: [PATCH 1/8] add perturbation to initial conditions in casadi solver --- pybamm/solvers/casadi_solver.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 87b80d4fde..111b1132e5 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -66,6 +66,10 @@ class CasadiSolver(pybamm.BaseSolver): return_solution_if_failed_early : bool, optional Whether to return a Solution object if the solver fails to reach the end of the simulation, but managed to take some successful steps. Default is False. + perturb_algebraic_initial_conditions : bool, optional + Whether to perturb algebraic initial conditions to avoid a singularity. This + can sometimes slow down the solver, but is kept True as default as it seems + to be more robust. """ def __init__( @@ -81,6 +85,7 @@ def __init__( extra_options_setup=None, extra_options_call=None, return_solution_if_failed_early=False, + perturb_algebraic_initial_conditions=True, ): super().__init__( "problem dependent", @@ -105,7 +110,7 @@ def __init__( self.extra_options_call = extra_options_call or {} self.extrap_tol = extrap_tol self.return_solution_if_failed_early = return_solution_if_failed_early - + self.perturb_algebraic_initial_conditions = perturb_algebraic_initial_conditions self.name = "CasADi solver with '{}' mode".format(mode) # Initialize @@ -666,6 +671,10 @@ def _run_integrator( y0_diff = y0[:len_rhs] y0_alg = y0[len_rhs:] + if self.perturb_algebraic_initial_conditions: + # Add a tiny perturbation to the algebraic initial conditions + # For some reason this helps with convergence + y0_alg = y0_alg * (1 + 1e-6 * np.random.rand(model.len_alg)) pybamm.logger.spam("Finished preliminary setup for integrator run") # Solve @@ -680,8 +689,9 @@ def _run_integrator( casadi_sol = integrator( x0=y0_diff, z0=y0_alg, p=inputs_with_tmin, **self.extra_options_call ) - except RuntimeError as e: + except RuntimeError as error: # If it doesn't work raise error + pybamm.logger.debug(f"Casadi integrator failed with error {error}") raise pybamm.SolverError(e.args[0]) pybamm.logger.debug("Finished casadi integrator") integration_time = timer.time() @@ -711,8 +721,9 @@ def _run_integrator( casadi_sol = integrator( x0=x, z0=z, p=inputs_with_tlims, **self.extra_options_call ) - except RuntimeError as e: + except RuntimeError as error: # If it doesn't work raise error + pybamm.logger.debug(f"Casadi integrator failed with error {error}") raise pybamm.SolverError(e.args[0]) integration_time = timer.time() x = casadi_sol["xf"] From 5332b7508f03b002b9f13b4689adb8d7da146e8a Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 11 Oct 2022 18:38:05 -0400 Subject: [PATCH 2/8] improve comment --- pybamm/solvers/casadi_solver.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 111b1132e5..5bc9536ab6 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -674,6 +674,8 @@ def _run_integrator( if self.perturb_algebraic_initial_conditions: # Add a tiny perturbation to the algebraic initial conditions # For some reason this helps with convergence + # The actual value of the initial conditions for the algebraic variables + # doesn't matter y0_alg = y0_alg * (1 + 1e-6 * np.random.rand(model.len_alg)) pybamm.logger.spam("Finished preliminary setup for integrator run") From 547df757748c330f9c745c52c489a1e6c75ff5d4 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 11 Oct 2022 18:39:02 -0400 Subject: [PATCH 3/8] flake8 --- pybamm/solvers/casadi_solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 5bc9536ab6..ad79a41a98 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -694,7 +694,7 @@ def _run_integrator( except RuntimeError as error: # If it doesn't work raise error pybamm.logger.debug(f"Casadi integrator failed with error {error}") - raise pybamm.SolverError(e.args[0]) + raise pybamm.SolverError(error.args[0]) pybamm.logger.debug("Finished casadi integrator") integration_time = timer.time() y_sol = casadi.vertcat(casadi_sol["xf"], casadi_sol["zf"]) @@ -726,7 +726,7 @@ def _run_integrator( except RuntimeError as error: # If it doesn't work raise error pybamm.logger.debug(f"Casadi integrator failed with error {error}") - raise pybamm.SolverError(e.args[0]) + raise pybamm.SolverError(error.args[0]) integration_time = timer.time() x = casadi_sol["xf"] z = casadi_sol["zf"] From 612ae75f969410371a1d5859c00ed3510f14496b Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 11 Oct 2022 22:47:14 -0400 Subject: [PATCH 4/8] fix tests --- pybamm/solvers/casadi_solver.py | 6 ++++-- .../test_experiments/test_simulation_with_experiment.py | 4 ++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index ad79a41a98..37a2c99482 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -663,20 +663,22 @@ def _run_integrator( integrator = self.integrators[model]["no grid"] len_rhs = model.concatenated_rhs.size + len_alg = model.concatenated_algebraic.size # Check y0 to see if it includes sensitivities if explicit_sensitivities: num_parameters = model.len_rhs_sens // model.len_rhs len_rhs = len_rhs * (num_parameters + 1) + len_alg = len_alg * (num_parameters + 1) y0_diff = y0[:len_rhs] y0_alg = y0[len_rhs:] - if self.perturb_algebraic_initial_conditions: + if self.perturb_algebraic_initial_conditions and len_alg > 0: # Add a tiny perturbation to the algebraic initial conditions # For some reason this helps with convergence # The actual value of the initial conditions for the algebraic variables # doesn't matter - y0_alg = y0_alg * (1 + 1e-6 * np.random.rand(model.len_alg)) + y0_alg = y0_alg * (1 + 1e-6 * casadi.DM(np.random.rand(len_alg))) pybamm.logger.spam("Finished preliminary setup for integrator run") # Solve diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index 8c24bac17b..a473842fee 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -476,11 +476,11 @@ def test_run_experiment_skip_steps(self): model, parameter_values=parameter_values, experiment=experiment2 ) sol2 = sim2.solve() - np.testing.assert_array_equal( + np.testing.assert_array_almost_equal( sol["Terminal voltage [V]"].data, sol2["Terminal voltage [V]"].data ) for idx1, idx2 in [(1, 0), (2, 1), (4, 2)]: - np.testing.assert_array_equal( + np.testing.assert_array_almost_equal( sol.cycles[0].steps[idx1]["Terminal voltage [V]"].data, sol2.cycles[0].steps[idx2]["Terminal voltage [V]"].data, ) From b21d89386c7497cc754940c14966af913a77a364 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Wed, 12 Oct 2022 10:09:39 -0400 Subject: [PATCH 5/8] no need to change tolerance in test --- CHANGELOG.md | 6 +++--- tests/integration/test_models/standard_model_tests.py | 5 ----- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a6c42188b..a3659b4742 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,14 +2,14 @@ ## Bug fixes -- For simulations with events that cause the simulation to stop early, the sensitivities could be evaluated incorrectly to zero ([#2331](https://github.com/pybamm-team/PyBaMM/pull/2337)) +- For simulations with events that cause the simulation to stop early, the sensitivities could be evaluated incorrectly to zero ([#2337](https://github.com/pybamm-team/PyBaMM/pull/2337)) ## Optimizations +- Added small perturbation to initial conditions for casadi solver. This seems to help the solver converge better in some cases ([#2356](https://github.com/pybamm-team/PyBaMM/pull/2356)) +- Added `ExplicitTimeIntegral` functionality to move variables which do not appear anywhere on the rhs to a new location, and to integrate those variables explicitly when `get` is called by the solution object. ([#2348](https://github.com/pybamm-team/PyBaMM/pull/2348)) - Added more rules for simplifying expressions ([#2211](https://github.com/pybamm-team/PyBaMM/pull/2211)) - - Sped up calculations of Electrode SOH variables for summary variables ([#2210](https://github.com/pybamm-team/PyBaMM/pull/2210)) -- Added `ExplicitTimeIntegral` functionality to move variables which do not appear anywhere on the rhs to a new location, and to integrate those variables explicitly when `get` is called by the solution object. ([#2348](https://github.com/pybamm-team/PyBaMM/pull/2348)) ## Breaking change diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py index b36a9ca16a..9d8040e516 100644 --- a/tests/integration/test_models/standard_model_tests.py +++ b/tests/integration/test_models/standard_model_tests.py @@ -72,11 +72,6 @@ def test_solving( self.solver.rtol = 1e-8 self.solver.atol = 1e-8 - # Somehow removing an equation makes the solver fail at - # the low tolerances - if isinstance(self.model, pybamm.lithium_ion.NewmanTobias): - self.solver.rtol = 1e-7 - Crate = abs( self.parameter_values["Current function [A]"] / self.parameter_values["Nominal cell capacity [A.h]"] From 3a7db9f3b381c796b41a3ca2a706918d4caca81e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 12 Oct 2022 14:10:52 +0000 Subject: [PATCH 6/8] style: pre-commit fixes --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index eaea0cf618..8027d4fa04 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,9 @@ [![black code style](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/ambv/black) + [![All Contributors](https://img.shields.io/badge/all_contributors-46-orange.svg)](#-contributors) + From a994f306345a8db0e84b94997f6d78b95e3d4178 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Wed, 12 Oct 2022 15:31:51 -0400 Subject: [PATCH 7/8] keep perturbation off for fast mode --- pybamm/solvers/casadi_solver.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 37a2c99482..7b4ef000bf 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -68,8 +68,8 @@ class CasadiSolver(pybamm.BaseSolver): the simulation, but managed to take some successful steps. Default is False. perturb_algebraic_initial_conditions : bool, optional Whether to perturb algebraic initial conditions to avoid a singularity. This - can sometimes slow down the solver, but is kept True as default as it seems - to be more robust. + can sometimes slow down the solver, but is kept True as default for "safe" mode + as it seems to be more robust (False by default for other modes). """ def __init__( @@ -85,7 +85,7 @@ def __init__( extra_options_setup=None, extra_options_call=None, return_solution_if_failed_early=False, - perturb_algebraic_initial_conditions=True, + perturb_algebraic_initial_conditions=None, ): super().__init__( "problem dependent", @@ -110,7 +110,18 @@ def __init__( self.extra_options_call = extra_options_call or {} self.extrap_tol = extrap_tol self.return_solution_if_failed_early = return_solution_if_failed_early - self.perturb_algebraic_initial_conditions = perturb_algebraic_initial_conditions + + # Decide whether to perturb algebraic initial conditions, True by default for + # "safe" mode, False by default for other modes + if perturb_algebraic_initial_conditions is None: + if mode == "safe": + self.perturb_algebraic_initial_conditions = True + else: + self.perturb_algebraic_initial_conditions = False + else: + self.perturb_algebraic_initial_conditions = ( + perturb_algebraic_initial_conditions + ) self.name = "CasADi solver with '{}' mode".format(mode) # Initialize From 647d3503bd0368b99d88d05f6d2d00369b87e513 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Wed, 12 Oct 2022 16:04:22 -0400 Subject: [PATCH 8/8] coverage --- tests/unit/test_solvers/test_casadi_solver.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py index 30923c40c2..6b9d073ab6 100644 --- a/tests/unit/test_solvers/test_casadi_solver.py +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -24,7 +24,12 @@ def test_model_solver(self): disc = pybamm.Discretisation() model_disc = disc.process_model(model, inplace=False) # Solve - solver = pybamm.CasadiSolver(mode="fast", rtol=1e-8, atol=1e-8) + solver = pybamm.CasadiSolver( + mode="fast", + rtol=1e-8, + atol=1e-8, + perturb_algebraic_initial_conditions=False, # added for coverage + ) t_eval = np.linspace(0, 1, 100) solution = solver.solve(model_disc, t_eval) np.testing.assert_array_equal(solution.t, t_eval)