Allow nonlinear constraints #353

[MarkBlyth]

Description

Enables nonlinear optimisation constraints. This turned out to be somewhat trivial with the newest PyBOP API; all it needed was a slightly different callback signature for the "trust-constr" optimiser.

Issue reference

Fixes #353

Review

Before you mark your PR as ready for review, please ensure that you've considered the following:

• Updated the CHANGELOG.md in reverse chronological order (newest at the top) with a concise description of the changes, including the PR number.
• Noted any breaking changes, including details on how it might impact existing functionality.

Type of change

• New Feature: A non-breaking change that adds new functionality.
• Optimization: A code change that improves performance.
• Examples: A change to existing or additional examples.
• Bug Fix: A non-breaking change that addresses an issue.
• Documentation: Updates to documentation or new documentation for new features.
• Refactoring: Non-functional changes that improve the codebase.
• Style: Non-functional changes related to code style (formatting, naming, etc).
• Testing: Additional tests to improve coverage or confirm functionality.
• Other: (Insert description of change)

Key checklist:

• No style issues: $ pre-commit run (or $ nox -s pre-commit) (see CONTRIBUTING.md for how to set this up to run automatically when committing locally, in just two lines of code)
• All unit tests pass: $ nox -s tests
• The documentation builds: $ nox -s doctest

You can run integration tests, unit tests, and doctests together at once, using $ nox -s quick.

Further checks:

• Code is well-commented, especially in complex or unclear areas.
• Added tests that prove my fix is effective or that my feature works.
• Checked that coverage remains or improves, and added tests if necessary to maintain or increase coverage.

Thank you for contributing to our project! Your efforts help us to deliver great software.

[NicolaCourtier]

Hi @MarkBlyth, great to hear that you could achieve the desired functionality with the current code structure!

I wonder if you could achieve the same result more efficiently (and for all optimisers) by extending the check_params function, because this function is run before the model is called rather than after. Did you try this?

[MarkBlyth]

I wonder if you could achieve the same result more efficiently (and for all optimisers) by extending the check_params function, because this function is run before the model is called rather than after. Did you try this?

Interesting idea, I hadn't thought of that! I've included an example now of using a subclassed Thevenin model, where the check_params does the constraining. Happy to integrate it into the core ECM file, if that's a route we want to go down.

I'd lean towards wanting to pass the parameter constraints into the optimiser, rather than the model. But given that this works easily and has no breaking changes, then perhaps it's the right way to go.

[NicolaCourtier]

Hi @MarkBlyth, thanks for the two examples! I've just merged #388 into develop so please check if it helps. Have you compared the performance of the two methods?

I like the extension of the SciPy callback function to wrap around base_callback as this would allow any other type of callback to be used as well, so we should add this feature.

For constraints though, I think the check_params method is more extensible and could be implemented similar to the SciPy method, i.e. I think we can pass a constraints function to be run within the check_params function, without defining a new class.

[MarkBlyth]

Have you compared the performance of the two methods?

trust-constr seems a bit faster and more robust.

The extra speed is possibly because it's a gradient-based method, which is suitable for this particular flavour of problem.

For robustness, I've found the gradient-free / metaheuristic optimisers can sometimes give up (for lack of better phrasing!) when I try to apply fairly tight bounds, whereas trust-constr manages. These tight bounds can be useful when trying to resolve fast timescales in GITT data. Maybe implementing the bounds at an optimiser level could help with that, e.g. changing sample-point distributions to ensure they always remain inside the bounds, so that the optimiser doesn't give up after seeing a big list of inf returns.

I think we can pass a constraints function to be run within the check_params function, without defining a new class.

Yep agreed, that sounds like the nicest way to do it. I'll get on that.

[MarkBlyth]

Finally got round to committing these changes. So far I've only set up the empirical models to take a user-defined parameter checker, but I'm sure the same could be done easily for physics-based models too

[NicolaCourtier]

Hi @MarkBlyth, looks like you're almost ready to merge but I noticed a small mistake in the develop version of the SciPyMinimize logging. Please could you consider fixing it alongside your additions?

The fix is to add a default value of self._cost0 = 1.0 to the SciPyMinimize.__init__ and line 197 to become:
self.log["cost"].append((cost if self.minimising else -cost)*self._cost0)
to rescale the cost value before logging. Otherwise I can create a separate issue and pull request...

[MarkBlyth]

@NicolaCourtier fix pushed, should be ready to merge now!

[BradyPlanden]

This is looking great, a few suggestions/requests. Happy to merge once they are resolved! Nice one @MarkBlyth

[BradyPlanden]

Thanks for the notebook! It looks great. Interestingly, it seems that the non-linear constraint removes the parameter bounds within the SciPy optimisation. In the parameter convergence plots, R1 is negative at the beginning (exceeds the lower bound previously defined) with an increasing cost. Any ideas as to why the bounds are not being applied to this optimisation?

[MarkBlyth]

I'm pretty sure the bounds are still being applied, but the default is to apply themwithout keep_feasible. PyBOP defines the bounds as a list (scipy_optimisers.py line 58)

           self._scipy_bounds = [
                (lower, upper)
                for lower, upper in zip(self.bounds["lower"], self.bounds["upper"])
            ]

but there's also the option to build a Bounds object instead

            self._scipy_bounds = Bounds(self.bounds["lower"], self.bounds["upper"], True)

The benefit of the Bounds object is that it offers a keep_feasible argument (that last True argument). This isn't an option when bounds are passed as a list, so keep_feasible will default to False for the standard PyBOP implementation, giving the negative resistances that we're seeing.

I'd lean towards keeping things as they are. The bounds are checked elsewhere in a way that's marginally more convenient with the current setup, and generally, the optimisers will perform better if they're not forced to stay within the feasible region. For the cases like this, where the infeasible parameters lead to simulation failure, the rest of the PyBOP infrastructure is able handles it properly to makes sure it doesn't lead to problems.

[BradyPlanden]

Thanks for the info on SciPy's keep_feasible argument, I wasn't aware of it. I've updated the SciPy bounds implementation to enable this functionality by default. The way I see it, if a user has requested bounds, they should be enforced, since we provide the functionality to not apply bounds on construction. While PyBOP has functionality to catch this in other areas, enforcing at the right level simplifies the complexity in the codebase. Please have a look at the change in the notebook and let me know if anything is wrong, the lower bound of the nonlinear constraint had to be reduced as the initial conditions were outside this value.

Overall, the performance appears to be improved with the keep_feasible option enabled for this problem.