Skip to content

Commit

Permalink
iRProp+ Optimiser (#636)
Browse files Browse the repository at this point in the history
* feat: adds iRProp+ optimiser

* removes redundant checks implemented in parent class, adds coverage

* fix: avoids potential multistart infeasible condition

* restores examples, updates to tell() comments

* tests: updates parameter range to align with default pybamm parameter object

* tests: updates priors for x0 within bound clipping for RandomSearch test

* docs: adds additional information to the iRProp+ docstring

* adds changelog entry
  • Loading branch information
BradyPlanden authored Feb 4, 2025
1 parent 9aba90f commit 5e19d0f
Show file tree
Hide file tree
Showing 8 changed files with 323 additions and 36 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
48 changes: 34 additions & 14 deletions examples/scripts/comparison_examples/irpropmin.py
Original file line number Diff line number Diff line change
@@ -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")
Expand Down
2 changes: 2 additions & 0 deletions pybop/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand All @@ -152,6 +153,7 @@
GradientDescent,
CMAES,
IRPropMin,
IRPropPlus,
NelderMead,
PSO,
SNES,
Expand Down
2 changes: 1 addition & 1 deletion pybop/optimisers/_adamw.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down
188 changes: 188 additions & 0 deletions pybop/optimisers/_irprop_plus.py
Original file line number Diff line number Diff line change
@@ -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
72 changes: 71 additions & 1 deletion pybop/optimisers/pints_optimisers.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
BasePintsOptimiser,
CuckooSearchImpl,
GradientDescentImpl,
IRPropPlusImpl,
RandomSearchImpl,
SimulatedAnnealingImpl,
)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down
Loading

0 comments on commit 5e19d0f

Please sign in to comment.