diff --git a/CHANGELOG.md b/CHANGELOG.md index 3e6cf962d2..6f2e22249f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Changed rootfinding algorithm to CasADi, scipy.optimize.root still accessible as an option ([#844](https://github.com/pybamm-team/PyBaMM/pull/844)) - Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824)) - Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819)) - Added options to export a solution to matlab or csv ([#811](https://github.com/pybamm-team/PyBaMM/pull/811)) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 1c031b8b8f..806015e900 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -23,7 +23,10 @@ class BaseSolver(object): atol : float, optional The absolute tolerance for the solver (default is 1e-6). root_method : str, optional - The method to use to find initial conditions (default is "lm") + The method to use to find initial conditions (default is "casadi"). If "casadi", + the solver uses casadi's Newton rootfinding algorithm to find initial + conditions. Otherwise, the solver uses 'scipy.optimize.root' with method + specified by 'root_method' (e.g. "lm", "hybr", ...) root_tol : float, optional The tolerance for the initial-condition solver (default is 1e-6). max_steps: int, optional @@ -36,7 +39,7 @@ def __init__( method=None, rtol=1e-6, atol=1e-6, - root_method="lm", + root_method="casadi", root_tol=1e-6, max_steps=1000, ): @@ -125,10 +128,11 @@ def set_up(self, model, inputs=None): "Cannot use ODE solver '{}' to solve DAE model".format(self.name) ) + if self.ode_solver is True: + self.root_method = None if ( - isinstance(self, pybamm.CasadiSolver) - and model.convert_to_format != "casadi" - ): + isinstance(self, pybamm.CasadiSolver) or self.root_method == "casadi" + ) and model.convert_to_format != "casadi": pybamm.logger.warning( f"Converting {model.name} to CasADi for solving with CasADi solver" ) @@ -268,6 +272,16 @@ def report(string): model.terminate_events_eval = terminate_events_eval model.discontinuity_events_eval = discontinuity_events_eval + # Save CasADi functions for the CasADi solver + # Note: when we pass to casadi the ode part of the problem must be in explicit + # form so we pre-multiply by the inverse of the mass matrix + if self.root_method == "casadi" or isinstance(self, pybamm.CasadiSolver): + mass_matrix_inv = casadi.MX(model.mass_matrix_inv.entries) + explicit_rhs = mass_matrix_inv @ rhs(t_casadi, y_casadi, u_casadi_stacked) + model.casadi_rhs = casadi.Function( + "rhs", [t_casadi, y_casadi, u_casadi_stacked], [explicit_rhs] + ) + model.casadi_algebraic = algebraic # Calculate consistent initial conditions for the algebraic equations if len(model.algebraic) > 0: all_states = pybamm.NumpyConcatenation( @@ -278,26 +292,13 @@ def report(string): model.residuals_eval = residuals_eval model.jacobian_eval = jacobian_eval y0_guess = y0.flatten() - model.y0 = self.calculate_consistent_state(model, 0, y0_guess) + model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs) else: # can use DAE solver to solve ODE model model.residuals_eval = Residuals(rhs, "residuals", model) model.jacobian_eval = jac_rhs model.y0 = y0.flatten() - # Save CasADi functions for the CasADi solver - # Note: when we pass to casadi the ode part of the problem must be in explicit - # form so we pre-multiply by the inverse of the mass matrix - if model.convert_to_format == "casadi" and isinstance( - self, pybamm.CasadiSolver - ): - mass_matrix_inv = casadi.MX(model.mass_matrix_inv.entries) - explicit_rhs = mass_matrix_inv @ rhs(t_casadi, y_casadi, u_casadi_stacked) - model.casadi_rhs = casadi.Function( - "rhs", [t_casadi, y_casadi, u_casadi_stacked], [explicit_rhs] - ) - model.casadi_algebraic = algebraic - pybamm.logger.info("Finish solver set-up") def set_inputs(self, model, ext_and_inputs): @@ -319,7 +320,7 @@ def set_inputs(self, model, ext_and_inputs): if model.jacobian_eval: model.jacobian_eval.set_inputs(ext_and_inputs) - def calculate_consistent_state(self, model, time=0, y0_guess=None): + def calculate_consistent_state(self, model, time=0, y0_guess=None, inputs=None): """ Calculate consistent state for the algebraic equations through root-finding @@ -328,6 +329,12 @@ def calculate_consistent_state(self, model, time=0, y0_guess=None): ---------- model : :class:`pybamm.BaseModel` The model for which to calculate initial conditions. + time : float + The time at which to calculate the states + y0_guess : :class:`np.array` + Guess for the rootfinding + inputs : dict, optional + Any input parameters to pass to the model when solving Returns ------- @@ -336,65 +343,101 @@ def calculate_consistent_state(self, model, time=0, y0_guess=None): of the algebraic equations) """ pybamm.logger.info("Start calculating consistent states") - rhs = model.rhs_eval - algebraic = model.algebraic_eval - jac = model.jac_algebraic_eval if y0_guess is None: y0_guess = model.concatenated_initial_conditions.flatten() # Split y0_guess into differential and algebraic - len_rhs = rhs(time, y0_guess).shape[0] + len_rhs = model.rhs_eval(time, y0_guess).shape[0] y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs]) + inputs = inputs or {} - def root_fun(y0_alg): - "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)" - y0 = np.concatenate([y0_diff, y0_alg]) - out = algebraic(time, y0) - pybamm.logger.debug( - "Evaluating algebraic equations at t={}, L2-norm is {}".format( - time * model.timescale, np.linalg.norm(out) - ) + # Solve using casadi or scipy + if self.root_method == "casadi": + # Set up + u_stacked = casadi.vertcat(*[x for x in inputs.values()]) + u = casadi.MX.sym("u", u_stacked.shape[0]) + y_alg = casadi.MX.sym("y_alg", y0_alg_guess.shape[0]) + y = casadi.vertcat(y0_diff, y_alg) + alg_root = model.casadi_algebraic(time, y, u) + # Solve + # set error_on_fail to False and just check the final output is small + # enough + roots = casadi.rootfinder( + "roots", + "newton", + dict(x=y_alg, p=u, g=alg_root), + {"abstol": self.root_tol}, ) - return out - - if jac: - if issparse(jac(0, y0_guess)): + try: + y0_alg = roots(y0_alg_guess, u_stacked).full().flatten() + success = True + message = None + # Check final output + fun = model.casadi_algebraic( + time, casadi.vertcat(y0_diff, y0_alg), u_stacked + ) + except RuntimeError as err: + success = False + message = err.args[0] + fun = None + else: + algebraic = model.algebraic_eval + jac = model.jac_algebraic_eval + + def root_fun(y0_alg): + "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)" + y0 = np.concatenate([y0_diff, y0_alg]) + out = algebraic(time, y0) + pybamm.logger.debug( + "Evaluating algebraic equations at t={}, L2-norm is {}".format( + time * model.timescale, np.linalg.norm(out) + ) + ) + return out - def jac_fn(y0_alg): - """ - Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) - """ - y0 = np.concatenate([y0_diff, y0_alg]) - return jac(0, y0)[:, len_rhs:].toarray() + if jac: + if issparse(jac(0, y0_guess)): - else: + def jac_fn(y0_alg): + """ + Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) + """ + y0 = np.concatenate([y0_diff, y0_alg]) + return jac(0, y0)[:, len_rhs:].toarray() - def jac_fn(y0_alg): - """ - Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) - """ - y0 = np.concatenate([y0_diff, y0_alg]) - return jac(0, y0)[:, len_rhs:] + else: - else: - jac_fn = None - # Find the values of y0_alg that are roots of the algebraic equations - sol = optimize.root( - root_fun, - y0_alg_guess, - jac=jac_fn, - method=self.root_method, - tol=self.root_tol, - ) - # Return full set of consistent initial conditions (y0_diff unchanged) - y0_consistent = np.concatenate([y0_diff, sol.x]) + def jac_fn(y0_alg): + """ + Evaluates jacobian using y0_diff (fixed) and y0_alg (varying) + """ + y0 = np.concatenate([y0_diff, y0_alg]) + return jac(0, y0)[:, len_rhs:] - if sol.success and np.all(sol.fun < self.root_tol * len(sol.x)): + else: + jac_fn = None + # Find the values of y0_alg that are roots of the algebraic equations + sol = optimize.root( + root_fun, + y0_alg_guess, + jac=jac_fn, + method=self.root_method, + tol=self.root_tol, + ) + # Set outputs + y0_alg = sol.x + success = sol.success + fun = sol.fun + message = sol.message + + if success and np.all(fun < self.root_tol * len(y0_alg)): + # Return full set of consistent initial conditions (y0_diff unchanged) + y0_consistent = np.concatenate([y0_diff, y0_alg]) pybamm.logger.info("Finish calculating consistent initial conditions") return y0_consistent - elif not sol.success: + elif not success: raise pybamm.SolverError( - "Could not find consistent initial conditions: {}".format(sol.message) + "Could not find consistent initial conditions: {}".format(message) ) else: raise pybamm.SolverError( @@ -402,7 +445,7 @@ def jac_fn(y0_alg): Could not find consistent initial conditions: solver terminated successfully, but maximum solution error ({}) above tolerance ({}) """.format( - np.max(sol.fun), self.root_tol * len(sol.x) + np.max(fun), self.root_tol * len(y0_alg) ) ) @@ -514,7 +557,7 @@ def solve(self, model, t_eval, external_variables=None, inputs=None): ) else: t_eval_dimensionless = np.insert( - t_eval_dimensionless, dindex, [dtime - eps, dtime + eps], + t_eval_dimensionless, dindex, [dtime - eps, dtime + eps] ) end_indices.append(len(t_eval_dimensionless)) @@ -548,7 +591,10 @@ def solve(self, model, t_eval, external_variables=None, inputs=None): last_state = solution.y[:, -1] if len(model.algebraic) > 0: model.y0 = self.calculate_consistent_state( - model, t_eval_dimensionless[end_index], last_state + model, + t_eval_dimensionless[end_index], + last_state, + ext_and_inputs, ) else: model.y0 = last_state diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index ea809fda7c..34fd740438 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -47,7 +47,7 @@ def __init__( mode="safe", rtol=1e-6, atol=1e-6, - root_method="lm", + root_method="casadi", root_tol=1e-6, max_step_decrease_count=5, **extra_options, diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py index 795f3361b8..be4a1932cc 100644 --- a/pybamm/solvers/idaklu_solver.py +++ b/pybamm/solvers/idaklu_solver.py @@ -36,7 +36,7 @@ class IDAKLUSolver(pybamm.BaseSolver): """ def __init__( - self, rtol=1e-6, atol=1e-6, root_method="lm", root_tol=1e-6, max_steps=1000 + self, rtol=1e-6, atol=1e-6, root_method="casadi", root_tol=1e-6, max_steps=1000 ): if idaklu_spec is None: diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py index 5f9999cbb7..74620cafed 100644 --- a/pybamm/solvers/scikits_dae_solver.py +++ b/pybamm/solvers/scikits_dae_solver.py @@ -40,7 +40,7 @@ def __init__( method="ida", rtol=1e-6, atol=1e-6, - root_method="lm", + root_method="casadi", root_tol=1e-6, max_steps=1000, ): diff --git a/setup.py b/setup.py index db141d4f4d..fa13bd203a 100644 --- a/setup.py +++ b/setup.py @@ -492,7 +492,7 @@ def load_version(): # List of dependencies install_requires=[ "numpy>=1.16", - "scipy>=1.0", + "scipy>=1.3", "pandas>=0.24", "anytree>=2.4.3", "autograd>=1.2", diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py index ce44197a34..7ab105982f 100644 --- a/tests/unit/test_solvers/test_base_solver.py +++ b/tests/unit/test_solvers/test_base_solver.py @@ -1,6 +1,7 @@ # # Tests for the Base Solver class # +import casadi import pybamm import numpy as np from scipy.sparse import csr_matrix @@ -49,9 +50,16 @@ def test_ode_solver_fail_with_dae(self): def test_find_consistent_initial_conditions(self): # Simple system: a single algebraic equation class ScalarModel: - concatenated_initial_conditions = np.array([[2]]) - jac_algebraic_eval = None - timescale = 1 + def __init__(self): + self.concatenated_initial_conditions = np.array([[2]]) + self.jac_algebraic_eval = None + self.timescale = 1 + t = casadi.MX.sym("t") + y = casadi.MX.sym("y") + u = casadi.MX.sym("u") + self.casadi_algebraic = casadi.Function( + "alg", [t, y, u], [self.algebraic_eval(t, y)] + ) def rhs_eval(self, t, y): return np.array([]) @@ -59,18 +67,30 @@ def rhs_eval(self, t, y): def algebraic_eval(self, t, y): return y + 2 - solver = pybamm.BaseSolver() + solver = pybamm.BaseSolver(root_method="lm") model = ScalarModel() init_cond = solver.calculate_consistent_state(model) np.testing.assert_array_equal(init_cond, -2) + # with casadi + solver_with_casadi = pybamm.BaseSolver(root_method="casadi", root_tol=1e-12) + model = ScalarModel() + init_cond = solver_with_casadi.calculate_consistent_state(model) + np.testing.assert_array_equal(init_cond, -2) # More complicated system vec = np.array([0.0, 1.0, 1.5, 2.0]) class VectorModel: - concatenated_initial_conditions = np.zeros_like(vec) - jac_algebraic_eval = None - timescale = 1 + def __init__(self): + self.concatenated_initial_conditions = np.zeros_like(vec) + self.jac_algebraic_eval = None + self.timescale = 1 + t = casadi.MX.sym("t") + y = casadi.MX.sym("y", vec.size) + u = casadi.MX.sym("u") + self.casadi_algebraic = casadi.Function( + "alg", [t, y, u], [self.algebraic_eval(t, y)] + ) def rhs_eval(self, t, y): return y[0:1] @@ -81,6 +101,9 @@ def algebraic_eval(self, t, y): model = VectorModel() init_cond = solver.calculate_consistent_state(model) np.testing.assert_array_almost_equal(init_cond, vec) + # with casadi + init_cond = solver_with_casadi.calculate_consistent_state(model) + np.testing.assert_array_almost_equal(init_cond, vec) # With jacobian def jac_dense(t, y): @@ -102,9 +125,16 @@ def jac_sparse(t, y): def test_fail_consistent_initial_conditions(self): class Model: - concatenated_initial_conditions = np.array([2]) - jac_algebraic_eval = None - timescale = 1 + def __init__(self): + self.concatenated_initial_conditions = np.array([2]) + self.jac_algebraic_eval = None + self.timescale = 1 + t = casadi.MX.sym("t") + y = casadi.MX.sym("y") + u = casadi.MX.sym("u") + self.casadi_algebraic = casadi.Function( + "alg", [t, y, u], [self.algebraic_eval(t, y)] + ) def rhs_eval(self, t, y): return np.array([]) @@ -120,12 +150,19 @@ def algebraic_eval(self, t, y): "Could not find consistent initial conditions: The iteration is not making", ): solver.calculate_consistent_state(Model()) - solver = pybamm.BaseSolver() + solver = pybamm.BaseSolver(root_method="lm") with self.assertRaisesRegex( pybamm.SolverError, "Could not find consistent initial conditions: solver terminated", ): solver.calculate_consistent_state(Model()) + # with casadi + solver = pybamm.BaseSolver(root_method="casadi") + with self.assertRaisesRegex( + pybamm.SolverError, + "Could not find consistent initial conditions: .../casadi", + ): + solver.calculate_consistent_state(Model()) def test_time_too_short(self): solver = pybamm.BaseSolver() diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index 9683b5cbc6..703673ae8e 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -28,7 +28,7 @@ def test_ida_roberts_klu(self): disc = pybamm.Discretisation() disc.process_model(model) - solver = pybamm.IDAKLUSolver() + solver = pybamm.IDAKLUSolver(root_method="lm") t_eval = np.linspace(0, 3, 100) solution = solver.solve(model, t_eval) @@ -58,7 +58,7 @@ def test_set_atol(self): mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) disc = pybamm.Discretisation(mesh, model.default_spatial_methods) disc.process_model(model) - solver = pybamm.IDAKLUSolver() + solver = pybamm.IDAKLUSolver(root_method="lm") variable_tols = {"Electrolyte concentration": 1e-3} solver.set_atol_by_variable(variable_tols, model) @@ -76,7 +76,7 @@ def test_failures(self): disc = pybamm.Discretisation() disc.process_model(model) - solver = pybamm.IDAKLUSolver() + solver = pybamm.IDAKLUSolver(root_method="lm") t_eval = np.linspace(0, 3, 100) with self.assertRaisesRegex(pybamm.SolverError, "KLU requires the Jacobian"): diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py index e3b39791a9..905507c692 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -194,7 +194,7 @@ def test_model_solver_dae_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -214,7 +214,7 @@ def test_model_solver_dae_bad_ics_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -238,7 +238,7 @@ def test_model_solver_dae_events_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 5, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_less(solution.y[0], 1.5) @@ -283,7 +283,7 @@ def nonsmooth_mult(t): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") # create two time series, one without a time point on the discontinuity, # and one with @@ -355,7 +355,7 @@ def jacobian(t, y): model.jacobian = jacobian # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -372,7 +372,7 @@ def test_solve_ode_model_with_dae_solver_python(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) @@ -421,7 +421,7 @@ def test_model_step_dae_python(self): disc = get_discretisation_for_testing() disc.process_model(model) - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") # Step once dt = 1 @@ -514,7 +514,10 @@ def test_model_solver_dae_inputs_events(self): disc.process_model(model) # Solve - solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) + if form == "python": + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm") + else: + solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8) t_eval = np.linspace(0, 5, 100) solution = solver.solve(model, t_eval, inputs={"rate 1": 0.1, "rate 2": 2}) np.testing.assert_array_less(solution.y[0], 1.5)