Differential method from Axen et al. (2022) for OCP with hysteresis

[chtamar]

Description

• Fixes Another method for OCP with hysteresis #4808

• Motivation: This method has been proposed by Axen et al. https://doi.org/10.1016/j.est.2022.103985, and like the Wycisk approach, it employs an ODE system to track the evolution of the ocp between the empirical lithiation and delithiation branches of the hysteresis. Although the Wycisk and Axen methods perform similarly, the proposed implementation does not need to compute the differentiation self.phase_param.U(sto_surf, T_bulk).diff(sto_surf), which makes it faster, and it also allows to choose separate decay rates for the lithiation and delithiation branches of the hysteresis. Important to note that more than one expression is proposed in https://doi.org/10.1016/j.est.2022.103985 for the delithiation branch, expression (9) in that paper is used here.

• Summary of the change:

  1. Employing the module wycisk_ocp as starting point, the new method is implemented in the module axen_ocp, which has been added in pybamm/models/submodels/interface/open_circuit_potential
  2. The file pybamm/models/submodels/interface/open_circuit_potential/__init__.py has been modified to make the axen_ocp module available in the namespace and import the class AxenOpenCircuitPotential
  3. The files pybamm/models/full_battery_models/base_battery_model.py and pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py have been modified to include the model's ocp option "Axen"
  4. The file pybamm/parameters/lithium_ion_parameters.py has been modified to add the new function U_hysteresis_branch, which includes the effect of temperature in the hysteresis branches, and to define the parameters K_lith and K_delith
  5. The file pybamm/CITATIONS.bib has been modified to add the relevant citation

• No extra dependencies are required

Type of change

Please add a line in the relevant section of CHANGELOG.md to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

• New feature (non-breaking change which adds functionality)
• Optimization (back-end change that speeds up the code)
• Bug fix (non-breaking change which fixes an issue)

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 tests pass: $ python -m pytest (or $ nox -s tests)
• The documentation builds: $ python -m pytest --doctest-plus src (or $ nox -s doctests)

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

Further checks:

• Code is commented, particularly in hard-to-understand areas
• Tests added that prove fix is effective or that feature works

[rtimms]

Thanks, @chtamar, the implementation looks great! Just a couple of comments, and could you also add some tests, please? Let me know if you need help with where to get started with the tests.

[chtamar]

Hi @rtimms, I added the requested tests, and everything seems to be working fine.

Initially, some tests failed with the error:

pybamm.expression_tree.exceptions.DomainError: children must have same or empty domains, not ['positive particle size'] and ['negative particle size']

To solve this problem, instead of defining ocp_bulk = delith_ref_x_av + H_x_av * h_x_av, the following definition is used, which was taken from the wycisk_ocp module:

delith_ref_s_av = pybamm.size_average(delith_ref_x_av)
H_s_av = pybamm.size_average(H_x_av)
h_s_av = pybamm.size_average(h_x_av)

ocp_bulk = delith_ref_s_av + H_s_av * h_s_av

This modification seems to solve the failed tests.

[rtimms]

Thanks, this looks great! As discussed in #4332, I would like to unify all the hysteresis options in PyBaMM. To help with that, I suggest we make it so that in all models, the hysteresis state variable goes between -1 and 1, with 0 corresponding to "equilibrium". Could you update it so that -1<=h<=1 and add a comment to the docs explaining the difference in implementation from the original paper?

[chtamar]

Hi @rtimms, the modification has been made and now the state variable h varies between -1 and 1, with the value 0 corresponding to the average OCP. A figure is attached below comparing the hysteresis state in Axen and Wycisk.

Regarding the comment explaining the different implementation from the original paper, should I add this as a jupyter notebook in PyBaMM\docs\source\examples\notebooks\models?

[rtimms]

Great, thanks! I would add a comment directly in the docstring of the new class and then update the existing Wycisk hysteresis notebook to compare all three models.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 98.70%. Comparing base (e39d72e) to head (c44cc41).
Report is 3 commits behind head on develop.

Additional details and impacted files

@@           Coverage Diff            @@
##           develop    #4816   +/-   ##
========================================
  Coverage    98.70%   98.70%           
========================================
  Files          303      304    +1     
  Lines        23353    23432   +79     
========================================
+ Hits         23050    23129   +79     
  Misses         303      303           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[chtamar]

Hi @rtimms, in the current PR, there are three different decay rates defined within the function pybamm.parameters.lithium_ion_parameters.ParticleLithiumIonParameters._set_parameters

• hysteresis_decay
• hysteresis_decay_lithiation
• hysteresis_decay_delithiation

Would it be better to define hysteresis_decay as an individual function directly under the class ParticleLithiumIonParameters?:

def hysteresis_decay(self, lithiation=None):
if lithiation is None:
lithiation = ""
else:
lithiation = lithiation + " "

return pybamm.Parameter(
f"{self.phase_prefactor}{Domain} particle {lithiation}hysteresis decay rate"
)

The module wycisk_ocp may also need a minor modification due to this change.

[rtimms]

Yes, sounds good! Anything that moves us towards unifying the hysteresis models is good. Doesn’t need to be for this PR, but we should make it so that the current sigmoid OCP also uses the state variable h to switch branches, then all the models are either an algebraic or ordinary differential equation for h (plus some specific implementation details).

[chtamar]

Hi @rtimms, I added the docstring to the model axen_ocp explaining the modifications made to the original method, so h changes between -1 and 1.

Also, the three parameters hysteresis_decay, hysteresis_decay_lithiation and hysteresis_decay_delithiation defined in the previous commit have been unified into a single function with the name hysteresis_decay. This led to a minor modification in the module wycisk_ocp, from K = self.phase_param.hysteresis_decay to K = self.phase_param.hysteresis_decay().

The only request that is still to be done is updating the existing Wycisk hysteresis notebook to compare the three models. This may take few days, to be sure it is right. I am wondering whether updating this notebook could be split into a separate issue?

[rtimms]

Thanks, this looks great, just needs a changelog entry. Please go ahead and make a new issue to update the notebook and we can do it in a separate PR.

[chtamar]

Hi @rtimms, an entry has been added to the changelog and the new issue #4851 has been created to update the existing Wycisk hysteresis notebook.

[rtimms]

Thanks for all your work on this @chtamar , looks great!

[rtimms]

@all-contributors please add @chtamar for code, docs and tests

[allcontributors]

@rtimms

I've put up a pull request to add @chtamar! 🎉

[chtamar]

Thanks for all your work on this @chtamar , looks great!

Thanks @rtimms, I really enjoyed working on this.

[rtimms]

@chtamar can you resolve the merge conflicts? Then this is good to merge

[chtamar]

@chtamar can you resolve the merge conflicts? Then this is good to merge

Hi @rtimms, my apologies for missing this. Conflict solved.