Repository: pybamm-team/PyBaMM (python)

Issue:
Another method for OCP with hysteresis
### Description

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. The rate at which the ocp moves back and forth between the hysteresis branches is scaled by the applied current and a decay rate.

### Motivation

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)`, and it also allows to choose separate decay rates for the lithiation and delithiation branches of the hysteresis. This option can be implemented for the Wycisk approach as well.

Together with the Wycisk approach, this method could be another alternative to choose from in pybamm.

A figure is attached comparing Wycisk and Axen for the silicon's ocp, sweeping over the decay rate (K) and temperature (T), using the _Chen2020_composite_ parametrization and an OCP entropic change for silicon of 2e-3 V/K, for testing purposes. Here the lithiation and delithiation decay rates for Axen take the same value.

![Image](https://github.com/user-attachments/assets/62fd504d-5c2e-42a8-9108-24c6c2a6de23)

### Possible Implementation

The method can be readily implemented using the module `wycisk_ocp` as starting point.

The method would need the following input parameters:
- The hysteresis branches for lithiation and delithition, `ocp_lith` and `ocp_delith`, respectively.
- The hysteresis decay rate for lithiation `K_lith`
- The hysteresis decay rate for delithiation `K_delith`
- The initial value for the state variable `h_init`

Main features of the method:
- As in the Wycisk approach, the state variable would be `h`, taking values between 0 and 1
- The rhs expression of the state variable would be scaled by the _volumetric interfacial current density_ `i_vol` and the decay rates `K_lith` and `K_delith`
- The magnitude of the OCP hysteresis would be tracked with the variable `H = lith_ref - delith_ref`, where `lith_ref = ocp_lith(sto, T)` and `delith_ref = ocp_delith(sto, T)`
- The resulting OCP is calculated as `ocp_surf = delith_ref + H * h`

The method would be implemented in the new file **_pybamm\models\submodels\interface\open_circuit_potential\axen_ocp.py_**

And the following files would be modified:
- **_pybamm\models\submodels\interface\open_circuit_potential\\\_\_init\_\_.py_**: to make the `axen_ocp` module available in the namespace and import the class `AxenOpenCircuitPotential`
- **_pybamm\models\full_battery_models\base_battery_model.py_** and **_pybamm\models\full_battery_models\lithium_ion\base_lithium_ion_model.py_**: to included the model's ocp option _"Axen"_
- **_pybamm\parameters\lithium_ion_parameters.py_**: 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`
- **_pybamm\CITATIONS.bib_**: to add the relevant citation

I could contribute with the implementation if you think this is worth including in pybamm.

### Additional context

_No response_

Pull request:
Differential method from Axen et al. (2022) for OCP with hysteresis
# Description

- Fixes #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](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] 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:

- [x] No style issues: `$ pre-commit run` (or `$ nox -s pre-commit`) (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) 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`)
- [x] 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


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c29a75cfb..8bc2ae96d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+- Added `axen_ocp` module within submodel `interface.open_circuit_potential` to handle an OCP with hysteresis. ([#4816](https://github.com/pybamm-team/PyBaMM/pull/4816))
 - Creates a 'calc_esoh' property in battery models ([#4825](https://github.com/pybamm-team/PyBaMM/pull/4825))
 - Added 'get_summary_variables' to return dictionary of computed summary variables ([#4824](https://github.com/pybamm-team/PyBaMM/pull/4824))
 - Added support for particle size distributions combined with particle mechanics. ([#4807](https://github.com/pybamm-team/PyBaMM/pull/4807))
diff --git a/docs/source/api/models/submodels/interface/open_circuit_potential/axen_ocp.rst b/docs/source/api/models/submodels/interface/open_circuit_potential/axen_ocp.rst
new file mode 100644
index 000000000..8e0e201b0
--- /dev/null
+++ b/docs/source/api/models/submodels/interface/open_circuit_potential/axen_ocp.rst
@@ -0,0 +1,5 @@
+Axen Open Circuit Potential
+=============================
+
+.. autoclass:: pybamm.open_circuit_potential.AxenOpenCircuitPotential
+    :members:
diff --git a/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst b/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst
index 3a20bccb1..84854567b 100644
--- a/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst
+++ b/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst
@@ -8,3 +8,4 @@ Open-circuit potential models
   single_ocp
   msmr_ocp
   wycisk_ocp
+  axen_ocp
diff --git a/src/pybamm/CITATIONS.bib b/src/pybamm/CITATIONS.bib
index f7ec22f09..d289ed5db 100644
--- a/src/pybamm/CITATIONS.bib
+++ b/src/pybamm/CITATIONS.bib
@@ -47,6 +47,19 @@
   doi = {10.1007/s12532-018-0139-4},
 }
 
+@article{Axen2022,
+title = {Evaluation of hysteresis expressions in a lumped voltage prediction model of a NiMH battery system in stationary storage applications},
+journal = {Journal of Energy Storage},
+volume = {48},
+pages = {103985},
+year = {2022},
+issn = {2352-152X},
+doi = {https://doi.org/10.1016/j.est.2022.103985},
+url = {https://www.sciencedirect.com/science/article/pii/S2352152X22000329},
+author = {Jenny B\"{o}rjesson Ax\'{e}n and Henrik Ekstr\"{o}m and Erika Widenkvist Zetterstr\"{o}m and G\"{o}ran Lindbergh},
+keywords = {Battery modelling, OCV hysteresis, NiMH, Voltage Response Model, Physical model, Application verification},
+}
+
 @article{Baker2018,
   title={Multi-species, multi-reaction model for porous intercalation electrodes: Part I. Model formulation and a perturbation solution for low-scan-rate, linear-sweep voltammetry of a spinel lithium manganese oxide electrode},
   author={Baker, Daniel R and Verbrugge, Mark W},
diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py
index 0678df8bd..956d44d00 100644
--- a/src/pybamm/models/full_battery_models/base_battery_model.py
+++ b/src/pybamm/models/full_battery_models/base_battery_model.py
@@ -102,9 +102,10 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 reactions.
             * "open-circuit potential" : str
                 Sets the model for the open circuit potential. Can be "single"
-                (default), "current sigmoid", "Wycisk", or "MSMR". If "MSMR" then the "particle"
-                option must also be "MSMR". A 2-tuple can be provided for different
-                behaviour in negative and positive electrodes.
+                (default), "current sigmoid", "Wycisk", "Axen", or "MSMR".
+                If "MSMR" then the "particle" option must also be "MSMR".
+                A 2-tuple can be provided for different behaviour in negative
+                and positive electrodes.
             * "operating mode" : str
                 Sets the operating mode for the model. This determines how the current
                 is set. Can be:
@@ -273,7 +274,13 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "stress and reaction-driven",
             ],
             "number of MSMR reactions": ["none"],
-            "open-circuit potential": ["single", "current sigmoid", "MSMR", "Wycisk"],
+            "open-circuit potential": [
+                "single",
+                "current sigmoid",
+                "MSMR",
+                "Wycisk",
+                "Axen",
+            ],
             "operating mode": [
                 "current",
                 "voltage",
diff --git a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
index 3c9595eb2..2e95c5c41 100644
--- a/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
+++ b/src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
@@ -273,6 +273,9 @@ class BaseModel(pybamm.BaseBatteryModel):
                 elif ocp_option == "Wycisk":
                     pybamm.citations.register("Wycisk2022")
                     ocp_model = ocp_submodels.WyciskOpenCircuitPotential
+                elif ocp_option == "Axen":
+                    pybamm.citations.register("Axen2022")
+                    ocp_model = ocp_submodels.AxenOpenCircuitPotential
                 elif ocp_option == "MSMR":
                     ocp_model = ocp_submodels.MSMROpenCircuitPotential
                 self.submodels[f"{domain} {phase} open-circuit potential"] = ocp_model(
diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/__init__.py b/src/pybamm/models/submodels/interface/open_circuit_potential/__init__.py
index 731327404..621e49d06 100644
--- a/src/pybamm/models/submodels/interface/open_circuit_potential/__init__.py
+++ b/src/pybamm/models/submodels/interface/open_circuit_potential/__init__.py
@@ -3,5 +3,6 @@ from .single_ocp import SingleOpenCircuitPotential
 from .current_sigmoid_ocp import CurrentSigmoidOpenCircuitPotential
 from .msmr_ocp import MSMROpenCircuitPotential
 from .wycisk_ocp import WyciskOpenCircuitPotential
+from .axen_ocp import AxenOpenCircuitPotential
 
-__all__ = ['base_ocp', 'current_sigmoid_ocp', 'msmr_ocp', 'single_ocp', 'wycisk_ocp']
+__all__ = ['base_ocp', 'current_sigmoid_ocp', 'msmr_ocp', 'single_ocp', 'wycisk_ocp', 'axen_ocp']
diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/axen_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/axen_ocp.py
new file mode 100644
index 000000000..9b2bf709f
--- /dev/null
+++ b/src/pybamm/models/submodels/interface/open_circuit_potential/axen_ocp.py
@@ -0,0 +1,150 @@
+#
+# from Axen 2022
+#
+import pybamm
+from . import BaseOpenCircuitPotential
+
+
+class AxenOpenCircuitPotential(BaseOpenCircuitPotential):
+    """
+    Class for open-circuit potential (OCP) with hysteresis based on the approach outlined by Axen et al. https://doi.org/10.1016/j.est.2022.103985.
+    This approach tracks the evolution of the hysteresis using an ODE system scaled by the volumetric interfacial current density and
+    a decay rate parameter.
+
+    The original method defines the open-circuit potential as OCP = OCP_delithiation + (OCP_lithiation - OCP_delithiation) * h,
+    where h is the hysteresis state. By allowing h to vary between 0 and 1, the OCP approaches OCP_lithiation when h tends to 1
+    and OCP_delithiation when h tends to 0. These interval limits are achieved by defining the ODE system as dh/dt = S_lith * (1 - h)
+    for lithiation and dh/dt = S_delith * h for delithiation, where S_lith and S_delith are rate coefficients.
+
+    The original implementation is modified to ensure h varies between -1 and 1, thus assuming an unified framework within pybamm.
+    This new variation interval is obtained with the ODE system dh/dt = S_lith * (1 + h) for lithiation and dh/dt = S_delith * (1 - h)
+    for delithiation. Now, OCP = OCP_delithiation + (OCP_lithiation - OCP_delithiation) * (1 - h) / 2, implying that the OCP approaches
+    OCP_lithiation when h tends to -1 and OCP_delithiation when h tends to 1.
+    """
+
+    def get_fundamental_variables(self):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        h = pybamm.Variable(
+            f"{Domain} electrode {phase_name}hysteresis state",
+            domains={
+                "primary": f"{domain} electrode",
+                "secondary": "current collector",
+            },
+        )
+        return {
+            f"{Domain} electrode {phase_name}hysteresis state": h,
+        }
+
+    def get_coupled_variables(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+
+        if self.reaction == "lithium-ion main":
+            T = variables[f"{Domain} electrode temperature [K]"]
+            h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+
+            # For "particle-size distribution" models, take distribution version
+            # of c_s_surf that depends on particle size.
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
+                sto_surf = variables[
+                    f"{Domain} {phase_name}particle surface stoichiometry distribution"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    sto_surf = sto_surf.orphans[0]
+                    T = T.orphans[0]
+                T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
+                h = pybamm.PrimaryBroadcast(h, [f"{domain} {phase_name}particle size"])
+            else:
+                sto_surf = variables[
+                    f"{Domain} {phase_name}particle surface stoichiometry"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    sto_surf = sto_surf.orphans[0]
+                    T = T.orphans[0]
+
+            variables[
+                f"{Domain} electrode {phase_name}hysteresis state distribution"
+            ] = h
+
+            # Bulk OCP is from the average SOC and temperature
+            lith_ref = self.phase_param.U(sto_surf, T, "lithiation")
+            delith_ref = self.phase_param.U(sto_surf, T, "delithiation")
+            H = lith_ref - delith_ref
+            variables[f"{Domain} electrode {phase_name}OCP hysteresis [V]"] = H
+
+            delith_ref_x_av = pybamm.x_average(delith_ref)
+            H_x_av = pybamm.x_average(H)
+            h_x_av = pybamm.x_average(h)
+            variables[f"X-averaged {domain} electrode {phase_name}hysteresis state"] = (
+                h_x_av
+            )
+
+            # check if psd
+            if domain_options["particle size"] == "distribution":
+                # should always be true
+                if f"{domain} particle size" in sto_surf.domains["primary"]:
+                    # check if MPM Model
+                    if "current collector" in sto_surf.domains["secondary"]:
+                        ocp_surf = delith_ref_x_av + H_x_av * (1 - h_x_av) / 2
+                    # must be DFN with PSD model
+                    elif (
+                        f"{domain} electrode" in sto_surf.domains["secondary"]
+                        or f"{domain} {phase_name}particle size"
+                        in sto_surf.domains["primary"]
+                    ):
+                        ocp_surf = delith_ref + H * (1 - h) / 2
+            # must not be a psd
+            else:
+                ocp_surf = delith_ref + H * (1 - h) / 2
+
+            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 * (1 - h_s_av) / 2
+
+            dUdT = self.phase_param.dUdT(sto_surf)
+
+        variables.update(self._get_standard_ocp_variables(ocp_surf, ocp_bulk, dUdT))
+        return variables
+
+    def set_rhs(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+
+        i_vol = variables[
+            f"{Domain} electrode {phase_name}volumetric interfacial current density [A.m-3]"
+        ]
+        c_max = self.phase_param.c_max
+        epsl = self.phase_param.epsilon_s
+
+        K_lith = self.phase_param.hysteresis_decay("lithiation")
+        K_delith = self.phase_param.hysteresis_decay("delithiation")
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+
+        S_lith = i_vol * K_lith / (self.param.F * c_max * epsl)
+        S_delith = i_vol * K_delith / (self.param.F * c_max * epsl)
+
+        i_vol_sign = pybamm.sign(i_vol)
+        signed_h = 1 - i_vol_sign * h
+        rate_coefficient = (
+            0.5 * (1 - i_vol_sign) * S_lith + 0.5 * (1 + i_vol_sign) * S_delith
+        )
+
+        dhdt = rate_coefficient * signed_h
+
+        self.rhs[h] = dhdt
+
+    def set_initial_conditions(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+        self.initial_conditions[h] = self.phase_param.h_init
diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py
index ae575589f..3e39a5740 100644
--- a/src/pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py
+++ b/src/pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py
@@ -160,7 +160,7 @@ class WyciskOpenCircuitPotential(BaseOpenCircuitPotential):
             f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
         ]
         dQdU = dQdU.orphans[0]
-        K = self.phase_param.hysteresis_decay
+        K = self.phase_param.hysteresis_decay()
         K_x = self.phase_param.hysteresis_switch
         h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
 
diff --git a/src/pybamm/parameters/lithium_ion_parameters.py b/src/pybamm/parameters/lithium_ion_parameters.py
index 118401953..6485047bf 100644
--- a/src/pybamm/parameters/lithium_ion_parameters.py
+++ b/src/pybamm/parameters/lithium_ion_parameters.py
@@ -454,9 +454,6 @@ class ParticleLithiumIonParameters(BaseParameters):
             self.epsilon_s * pybamm.r_average(self.c_init)
         )
         # if self.options['open-circuit potential'] == 'Plett':
-        self.hysteresis_decay = pybamm.Parameter(
-            f"{pref}{Domain} particle hysteresis decay rate"
-        )
         self.hysteresis_switch = pybamm.Parameter(
             f"{pref}{Domain} particle hysteresis switching factor"
         )
@@ -513,6 +510,21 @@ class ParticleLithiumIonParameters(BaseParameters):
         self.b_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant b")
         self.m_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant m")
 
+    def hysteresis_decay(self, lithiation=None):
+        """
+        Rate at which the open-circuit potential approaches the lithiation
+        or delithiation branch when it exhibits hysteresis.
+        """
+        Domain = self.domain.capitalize()
+        if lithiation is None:
+            lithiation = ""
+        else:
+            lithiation = lithiation + " "
+
+        return pybamm.Parameter(
+            f"{self.phase_prefactor}{Domain} particle {lithiation}hysteresis decay rate"
+        )
+
     def k_cr(self, T):
         """
         Cracking rate for the electrode;


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
index 1a093c68e..8764ec510 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
@@ -90,6 +90,30 @@ class TestMPM:
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all(skip_output_tests=True)
 
+    def test_axen_ocp(self):
+        options = {"open-circuit potential": ("Axen", "single")}
+        model = pybamm.lithium_ion.MPM(options)
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values = pybamm.get_size_distribution_parameters(parameter_values)
+        parameter_values.update(
+            {
+                "Negative electrode lithiation OCP [V]": lambda sto: parameter_values[
+                    "Negative electrode OCP [V]"
+                ](sto)
+                - 0.1,
+                "Negative electrode delithiation OCP [V]": lambda sto: parameter_values[
+                    "Negative electrode OCP [V]"
+                ](sto)
+                + 0.1,
+                "Negative particle lithiation hysteresis decay rate": 10,
+                "Negative particle delithiation hysteresis decay rate": 10,
+                "Initial hysteresis state in negative electrode": 0.0,
+            },
+            check_already_exists=False,
+        )
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all(skip_output_tests=True)
+
     def test_voltage_control(self):
         options = {"operating mode": "voltage"}
         model = pybamm.lithium_ion.MPM(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index dc056523a..2950f1e93 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -34,7 +34,7 @@ PRINT_OPTIONS_OUTPUT = """\
 'lithium plating porosity change': 'false' (possible: ['false', 'true'])
 'loss of active material': 'stress-driven' (possible: ['none', 'stress-driven', 'reaction-driven', 'current-driven', 'stress and reaction-driven'])
 'number of MSMR reactions': 'none' (possible: ['none'])
-'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR', 'Wycisk'])
+'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR', 'Wycisk', 'Axen'])
 'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'differential power', 'explicit power', 'resistance', 'differential resistance', 'explicit resistance', 'CCCV'])
 'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'uniform profile', 'quadratic profile', 'quartic profile', 'MSMR'])
 'particle mechanics': 'swelling only' (possible: ['none', 'swelling only', 'swelling and cracking'])
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index d7d061f2b..010e59aa6 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -503,6 +503,10 @@ class BaseUnitTestLithiumIon:
         options = {"open-circuit potential": "Wycisk"}
         self.check_well_posedness(options)
 
+    def test_well_posed_axen_ocp(self):
+        options = {"open-circuit potential": "Axen"}
+        self.check_well_posedness(options)
+
     def test_well_posed_msmr(self):
         options = {
             "open-circuit potential": "MSMR",
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 56f4c76f4..e5f5b11c7 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -52,6 +52,20 @@ class TestDFN(BaseUnitTestLithiumIon):
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_axen_ocp_with_psd(self):
+        options = {
+            "open-circuit potential": "Axen",
+            "particle size": "distribution",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_axen_ocp_with_composite(self):
+        options = {
+            "open-circuit potential": (("Axen", "single"), "single"),
+            "particle phases": ("2", "1"),
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_external_circuit_explicit_power(self):
         options = {"operating mode": "explicit power"}
         self.check_well_posedness(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
index c548aba49..d28eb2374 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
@@ -103,6 +103,13 @@ class TestMPM:
         model = pybamm.lithium_ion.MPM(options)
         model.check_well_posedness()
 
+    def test_axen_ocp(self):
+        options = {
+            "open-circuit potential": "Axen",
+        }
+        model = pybamm.lithium_ion.MPM(options)
+        model.check_well_posedness()
+
     def test_mpm_with_lithium_plating(self):
         options = {
             "lithium plating": "irreversible",