diff --git a/CHANGELOG.md b/CHANGELOG.md index f58d0485..6acdf807 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#636](https://github.com/pybop-team/PyBOP/pull/636) - Adds `pybop.IRPropPlus` optimiser with corresponding tests. - [#635](https://github.com/pybop-team/PyBOP/pull/635) - Adds support for multi-proposal evaluation of list-like objects to `BaseCost` classes. - [#635](https://github.com/pybop-team/PyBOP/pull/635) - Adds global parameter sensitivity analysis with method `BaseCost.sensitivity_analysis`. This is computation is added to `OptimisationResult` if optimiser arg `compute_sensitivities` is `True`. An additional arg is added to select the number of samples for analysis: `n_sensitivity_samples`. - [#630] (https://github.com/pybop-team/PyBOP/pull/632) - Fisher Information Matrix added to `BaseLikelihood` class. diff --git a/examples/scripts/comparison_examples/irpropmin.py b/examples/scripts/comparison_examples/irpropmin.py index 403643e4..230fd693 100644 --- a/examples/scripts/comparison_examples/irpropmin.py +++ b/examples/scripts/comparison_examples/irpropmin.py @@ -1,46 +1,66 @@ import numpy as np -import pybamm import pybop # Define model and use high-performant solver for sensitivities -solver = pybamm.IDAKLUSolver() parameter_set = pybop.ParameterSet("Chen2020") -model = pybop.lithium_ion.SPM(parameter_set=parameter_set, solver=solver) +model = pybop.lithium_ion.SPM(parameter_set=parameter_set) -# Fitting parameters +# Construct the fitting parameters +# with initial values sampled from a different distribution parameters = pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.02), + prior=pybop.Uniform(0.3, 0.9), + bounds=[0.3, 0.8], + initial_value=pybop.Uniform(0.4, 0.75).rvs()[0], + true_value=parameter_set["Negative electrode active material volume fraction"], ), pybop.Parameter( "Positive electrode active material volume fraction", - prior=pybop.Gaussian(0.48, 0.02), + prior=pybop.Uniform(0.3, 0.9), + initial_value=pybop.Uniform(0.4, 0.75).rvs()[0], + true_value=parameter_set["Positive electrode active material volume fraction"], + # no bounds ), ) # Generate data -sigma = 0.001 -t_eval = np.arange(0, 900, 3) -values = model.predict(t_eval=t_eval) -corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval)) +sigma = 0.002 +experiment = pybop.Experiment( + [ + ( + "Discharge at 0.5C for 12 minutes (10 second period)", + "Charge at 0.5C for 12 minutes (10 second period)", + ) + ] +) +values = model.predict(initial_state={"Initial SoC": 0.4}, experiment=experiment) + + +def noise(sigma): + return np.random.normal(0, sigma, len(values["Voltage [V]"].data)) + # Form dataset dataset = pybop.Dataset( { - "Time [s]": t_eval, + "Time [s]": values["Time [s]"].data, "Current function [A]": values["Current [A]"].data, - "Voltage [V]": corrupt_values, + "Voltage [V]": values["Voltage [V]"].data + noise(sigma), } ) # Generate problem, cost function, and optimisation class problem = pybop.FittingProblem(model, parameters, dataset) -cost = pybop.Minkowski(problem, p=2) -optim = pybop.IRPropMin(cost, max_iterations=100, sigma0=0.1) +likelihood = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=sigma * 1.05) +posterior = pybop.LogPosterior(likelihood) +optim = pybop.IRPropMin( + posterior, max_iterations=125, max_unchanged_iterations=60, sigma0=0.01 +) results = optim.run() +print(parameters.true_value()) # Plot the timeseries output pybop.plot.quick(problem, problem_inputs=results.x, title="Optimised Comparison") diff --git a/pybop/__init__.py b/pybop/__init__.py index 8fff7052..2b60cce7 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -141,6 +141,7 @@ from .optimisers._adamw import AdamWImpl from .optimisers._gradient_descent import GradientDescentImpl from .optimisers._simulated_annealing import SimulatedAnnealingImpl +from .optimisers._irprop_plus import IRPropPlusImpl from .optimisers.base_optimiser import BaseOptimiser, OptimisationResult from .optimisers.base_pints_optimiser import BasePintsOptimiser from .optimisers.scipy_optimisers import ( @@ -152,6 +153,7 @@ GradientDescent, CMAES, IRPropMin, + IRPropPlus, NelderMead, PSO, SNES, diff --git a/pybop/optimisers/_adamw.py b/pybop/optimisers/_adamw.py index 678a57f2..b01d2dea 100644 --- a/pybop/optimisers/_adamw.py +++ b/pybop/optimisers/_adamw.py @@ -1,5 +1,5 @@ # -# Extends the Pints' Adam Class with a Weight Decay addition +# Initially based off Pints' Adam class, adds weight decay # import warnings diff --git a/pybop/optimisers/_irprop_plus.py b/pybop/optimisers/_irprop_plus.py new file mode 100644 index 00000000..a2ff4692 --- /dev/null +++ b/pybop/optimisers/_irprop_plus.py @@ -0,0 +1,188 @@ +# +# Initially based of Pints' IRProp- class. +# + +import numpy as np +from pints import Optimiser as PintsOptimiser +from pints import RectangularBoundaries + + +class IRPropPlusImpl(PintsOptimiser): + """ + The iRprop+ algorithm adjusts step sizes based on the sign of the gradient + in each dimension, increasing step sizes when the sign remains consistent + and decreasing when the sign changes. This implementation includes + weight (parameter) backtracking that reverts the previous step in the + event of a gradient sign changing. + + References: + - [1] Igel and Hüsken (2003): Empirical Evaluation of Improved Rprop Algorithms. + - [2] Riedmiller and Braun (1993): A Direct Adaptive Method for Faster Backpropagation. + - [3] Igel and Husk (2003): Improving the Rprop Learning Algorithm. + + Parameters + ---------- + x0 : array-like + Initial starting point for the optimisation. + sigma0 : float or array-like, optional + Initial step size(s). If a scalar is provided, it is applied to all dimensions. + Default is 0.05. + boundaries : pints.Boundaries, optional + Boundary constraints for the optimisation. If None, no boundaries are applied. + + Attributes + ---------- + eta_min : float + Factor by which the step size is reduced when the gradient sign changes. + Default is 0.5. + eta_max : float + Factor by which the step size is increased when the gradient sign remains consistent. + Default is 1.2. + step_min : float + Minimum allowable step size. Default is 1e-3 * min(sigma0). + step_max : float or None + Maximum allowable step size. Default is None (unlimited). + """ + + def __init__(self, x0, sigma0=0.05, boundaries=None): + super().__init__(x0, sigma0, boundaries) + + # Set hypers + self.eta_min = 0.5 + self.eta_max = 1.2 + self.step_min = 1e-4 * np.min(self._sigma0) + self.step_max = None + + # Store the previous update for backtracking + self._update_previous = np.zeros_like(x0, dtype=float) + + # Internal states + self._x_current = np.array(x0, dtype=float) + self._gradient_previous = None + self._step_sizes = np.copy(self._sigma0) + self._f_current = np.inf + self._x_best = np.array(x0, dtype=float) + self._f_best = np.inf + self._running = False + self._ready_for_tell = False + + # Boundaries + self._boundaries = boundaries + if isinstance(boundaries, RectangularBoundaries): + self._lower = boundaries.lower() + self._upper = boundaries.upper() + else: + self._lower = self._upper = None + + # Set proposals + self._proposed = np.copy(self._x_current) + self._proposed.setflags(write=False) + + def ask(self): + """ + Proposes the next point for evaluation. + + Returns + ------- + list + A list containing the proposed point. + """ + if not self._running: + if self.step_min is not None and self.step_max is not None: + if self.step_min >= self.step_max: + raise ValueError( + f"Minimum step size ({self.step_min}) must be less than " + f"maximum step size ({self.step_max})." + ) + self._running = True + + self._ready_for_tell = True + return [self._proposed] + + def tell(self, reply): + """ + Updates the optimiser with the function value and gradient at the proposed point. + + Parameters + ---------- + reply : list + A list containing a tuple of (function value, gradient) at the proposed point. + + Raises + ------ + RuntimeError + If `ask()` was not called before `tell()`. + """ + if not self._ready_for_tell: + raise RuntimeError("ask() must be called before tell().") + + self._ready_for_tell = False + f_new, gradient_new = reply[0] + + # Setup for first iteration + if self._gradient_previous is None: + self._gradient_previous = gradient_new + self._f_current = f_new + return + + # Compute element-wise gradient product + grad_product = gradient_new * self._gradient_previous + + # Update step sizes, and bound them + self._step_sizes[grad_product > 0] *= self.eta_max + self._step_sizes[grad_product < 0] *= self.eta_min + self._step_sizes = np.clip(self._step_sizes, self.step_min, self.step_max) + + # Handle weight backtracking, + # Reverting last update where gradient sign changed + gradient_new[grad_product < 0] = 0 + self._x_current[grad_product < 0] -= self._update_previous[grad_product < 0] + + # Update the current position + self._x_current = np.copy(self._proposed) + self._f_current = f_new + self._gradient_previous = gradient_new + + # Compute the new update and store for back-tracking + update = -self._step_sizes * np.sign(gradient_new) + self._update_previous = update + + # Step in the direction of the negative gradient + proposed = self._x_current + update + + # Boundaries + if self._lower is not None: + # Rectangular boundaries + while np.any(proposed < self._lower) or np.any(proposed >= self._upper): + mask = np.logical_or(proposed < self._lower, proposed >= self._upper) + self._step_sizes[mask] *= self.eta_min + proposed = self._x_current - self._step_sizes * np.sign(gradient_new) + + # Update proposed attribute + self._proposed = proposed + self._proposed.setflags(write=False) + + # Update best solution + if f_new < self._f_best: + self._f_best = f_new + self._x_best = self._x_current + + def running(self): + """Returns the state of the optimiser""" + return self._running + + def f_best(self): + """Returns the best function value found so far.""" + return self._f_best + + def x_best(self): + """Returns the best position found so far.""" + return self._x_best + + def name(self): + """Returns the name of the optimiser.""" + return "iRprop+" + + def needs_sensitivities(self): + """Indicates that this optimiser requires gradient information.""" + return True diff --git a/pybop/optimisers/pints_optimisers.py b/pybop/optimisers/pints_optimisers.py index 39841145..476cc543 100644 --- a/pybop/optimisers/pints_optimisers.py +++ b/pybop/optimisers/pints_optimisers.py @@ -10,6 +10,7 @@ BasePintsOptimiser, CuckooSearchImpl, GradientDescentImpl, + IRPropPlusImpl, RandomSearchImpl, SimulatedAnnealingImpl, ) @@ -164,7 +165,7 @@ class IRPropMin(BasePintsOptimiser): Implements the iRpropMin optimisation algorithm. This class inherits from the PINTS IRPropMin class, which is an optimiser that - uses resilient backpropagation with weight-backtracking. It is designed to handle + uses resilient backpropagation without weight-backtracking. It is designed to handle problems with large plateaus, noisy gradients, and local minima. Parameters @@ -229,6 +230,75 @@ def __init__( ) +class IRPropPlus(BasePintsOptimiser): + """ + Implements the iRpropPlus optimisation algorithm. + + This class implements the improved resilient backpropagation with weight-backtracking. + It is designed to handle problems with large plateaus, noisy gradients, and local minima. + + Parameters + ---------- + cost : callable + The cost function to be minimized. + max_iterations : int, optional + Maximum number of iterations for the optimisation. + min_iterations : int, optional (default=2) + Minimum number of iterations before termination. + max_unchanged_iterations : int, optional (default=15) + Maximum number of iterations without improvement before termination. + multistart : int, optional (default=1) + Number of optimiser restarts from randomly sample position. These positions + are sampled from the priors. + parallel : bool, optional (default=False) + Whether to run the optimisation in parallel. + **optimiser_kwargs : optional + Valid PINTS option keys and their values, for example: + x0 : array_like + Initial position from which optimisation will start. + sigma0 : float + Initial step size or standard deviation depending on the optimiser. + bounds : dict + A dictionary with 'lower' and 'upper' keys containing arrays for lower and + upper bounds on the parameters. + use_f_guessed : bool + Whether to return the guessed function values. + absolute_tolerance : float + Absolute tolerance for convergence checking. + relative_tolerance : float + Relative tolerance for convergence checking. + max_evaluations : int + Maximum number of function evaluations. + threshold : float + Threshold value for early termination. + + See Also + -------- + pints.IRPropMin : The PINTS implementation this class is based on. + """ + + def __init__( + self, + cost, + max_iterations: int = None, + min_iterations: int = 2, + max_unchanged_iterations: int = 15, + multistart: int = 1, + parallel: bool = False, + **optimiser_kwargs, + ): + super().__init__( + cost, + IRPropPlusImpl, + max_iterations, + min_iterations, + max_unchanged_iterations, + multistart, + parallel, + **optimiser_kwargs, + ) + + class PSO(BasePintsOptimiser): """ Implements a particle swarm optimisation (PSO) algorithm. diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index ae29fad2..25f64a0a 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -79,6 +79,7 @@ def noise(self, sigma, values): pybop.CuckooSearch, pybop.NelderMead, pybop.IRPropMin, + pybop.IRPropPlus, pybop.AdamW, pybop.CMAES, pybop.SNES, diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index 43053af2..91069cdd 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -31,8 +31,8 @@ def dataset(self): def one_parameter(self): return pybop.Parameter( "Positive electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.02), - bounds=[0.58, 0.62], + prior=pybop.Gaussian(0.5, 0.02), + bounds=[0.48, 0.52], ) @pytest.fixture @@ -40,13 +40,13 @@ def two_parameters(self): return pybop.Parameters( pybop.Parameter( "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.2), + prior=pybop.Gaussian(0.6, 0.02), bounds=[0.58, 0.62], ), pybop.Parameter( "Positive electrode active material volume fraction", - prior=pybop.Gaussian(0.55, 0.05), - bounds=[0.53, 0.57], + prior=pybop.Gaussian(0.5, 0.05), + bounds=[0.48, 0.52], ), ) @@ -89,6 +89,7 @@ def two_param_cost(self, model, two_parameters, dataset): (pybop.XNES, "Exponential Natural Evolution Strategy (xNES)", False), (pybop.PSO, "Particle Swarm Optimisation (PSO)", False), (pybop.IRPropMin, "iRprop-", True), + (pybop.IRPropPlus, "iRprop+", True), (pybop.NelderMead, "Nelder-Mead", False), (pybop.RandomSearch, "Random Search", False), (pybop.SimulatedAnnealing, "Simulated Annealing", False), @@ -137,6 +138,7 @@ def test_no_optimisation_parameters(self, model, dataset): pybop.XNES, pybop.PSO, pybop.IRPropMin, + pybop.IRPropPlus, pybop.NelderMead, pybop.CuckooSearch, pybop.RandomSearch, @@ -189,7 +191,7 @@ def check_multistart(optim, n_iters, multistarts): check_incorrect_update(optim) # Test multistart - multistart_optim = optimiser(cost, multistart=2, max_iterations=6) + multistart_optim = optimiser(cost, max_iterations=6, multistart=2) check_multistart(multistart_optim, 6, 2) if optimiser in [pybop.GradientDescent, pybop.AdamW, pybop.NelderMead]: @@ -201,7 +203,7 @@ def check_multistart(optim, n_iters, multistarts): optim, {"upper": [np.inf], "lower": [0.57]}, should_raise=True ) else: - bounds = {"upper": [0.63], "lower": [0.57]} + bounds = {"upper": [0.53], "lower": [0.47]} check_bounds_handling(optim, cost_bounds) optim = optimiser(cost=cost, bounds=bounds) assert optim.bounds == bounds @@ -282,6 +284,7 @@ def check_multistart(optim, n_iters, multistarts): if optimiser in [ pybop.AdamW, + pybop.IRPropPlus, pybop.CuckooSearch, pybop.GradientDescent, pybop.RandomSearch, @@ -370,7 +373,7 @@ def check_multistart(optim, n_iters, multistarts): else: x0 = cost.parameters.initial_value() assert optim.x0 == x0 - x0_new = np.array([0.6]) + x0_new = np.array([0.5]) optim = optimiser(cost=cost, x0=x0_new) assert optim.x0 == x0_new assert optim.x0 != x0 @@ -382,21 +385,17 @@ def test_cuckoo_no_bounds(self, cost): def test_randomsearch_bounds(self, two_param_cost): # Test clip_candidates with bound - bounds = {"upper": [0.62, 0.57], "lower": [0.58, 0.53]} - optimiser = pybop.RandomSearch( - cost=two_param_cost, bounds=bounds, max_iterations=1 - ) - candidates = np.array([[0.57, 0.52], [0.63, 0.58]]) - clipped_candidates = optimiser.optimiser.clip_candidates(candidates) - expected_clipped = np.array([[0.58, 0.53], [0.62, 0.57]]) + bounds = {"upper": [0.62, 0.54], "lower": [0.58, 0.46]} + optim = pybop.RandomSearch(cost=two_param_cost, bounds=bounds, max_iterations=1) + candidates = np.array([[0.57, 0.55], [0.63, 0.44]]) + clipped_candidates = optim.optimiser.clip_candidates(candidates) + expected_clipped = np.array([[0.58, 0.54], [0.62, 0.46]]) assert np.allclose(clipped_candidates, expected_clipped) # Test clip_candidates without bound - optimiser = pybop.RandomSearch( - cost=two_param_cost, bounds=None, max_iterations=1 - ) + optim = pybop.RandomSearch(cost=two_param_cost, bounds=None, max_iterations=1) candidates = np.array([[0.57, 0.52], [0.63, 0.58]]) - clipped_candidates = optimiser.optimiser.clip_candidates(candidates) + clipped_candidates = optim.optimiser.clip_candidates(candidates) assert np.allclose(clipped_candidates, candidates) def test_randomsearch_ask_without_bounds(self, two_param_cost): @@ -422,6 +421,12 @@ def test_adamw_impl_bounds(self): with pytest.warns(UserWarning, match="Boundaries ignored by AdamW"): pybop.AdamWImpl(x0=[0.1], boundaries=[0.0, 0.2]) + def test_irprop_plus_impl_incorrect_steps(self): + with pytest.raises(ValueError, match="Minimum step size"): + optim = pybop.IRPropPlusImpl(x0=[0.1]) + optim.step_max = 1e-8 + optim.ask() + def test_scipy_minimize_with_jac(self, cost): # Check a method that uses gradient information optim = pybop.SciPyMinimize( @@ -732,7 +737,7 @@ def test_optimisation_results(self, cost): # Test list-like functionality with "best" properties optim = pybop.Optimisation( cost=cost, - x0=[0.6, 0.5], + x0=[0.5], n_iterations=1, multistart=3, compute_sensitivities=True,