Skip to content

Issue 4884 unify hysteresis #4893

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 6 commits into
base: develop
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view

This file was deleted.

356 changes: 356 additions & 0 deletions docs/source/examples/notebooks/models/hysteresis-state-models.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,356 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Hysteresis State models"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.2.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m25.0.1\u001b[0m\n",
"\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n",
"Note: you may need to restart the kernel to use updated packages.\n"
]
}
],
"source": [
"%pip install \"pybamm[plot,cite]\" -q # install PyBaMM if it is not installed\n",
"import pybamm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Model Equations\n",
"\n",
"Herein, we outline the equations for the \"Curret Sigmoid\" model, as described in Ai et al. (2022), the Differential Capacity Hysteresis State open-circuit potential model, as described in Wycisk (2022), and the Hysteresis State open-circuit potential model, as described in Axen (2022).\n",

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo: "Curret" -> "Current"

"\n",
"In all of the models, the open-circuit potential is given by\n",
"\n",
"$$U = \\frac{1+h}{2} U_{delith} + \\frac{1-h}{2} U_{lith},$$\n",

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Prefer throughout to use Roman (non-italic) text for textual subscripts.

"\n",
"where $h$ is a hysteresis state variable, $U_{delith}$ is the delithiation open-circuit potential, and $U_{lith}$ is the lithiation open-circuit potential. The hysteresis state variable $h$ is defined between -1 and 1, where -1 corresponds to the delithiation branch and 1 corresponds to the lithiation branch.\n",

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definition of $h$ is back-to-front in text vs equation.

"\n",
"### Current Sigmoid (Ai)\n",
"This approach uses a sigmoid function to switch between the delithiation and lithiation branches of the open-circuit potential, depending on the sign of the current. The original paper gives the following expression for the open-circuit potential:\n",
"\n",
"$$ U = \\text{sigmoid}\\left(-\\frac{KI}{Q_{cell}}\\right) U_{delith} + \\text{sigmoid}\\left(\\frac{KI}{Q_{cell}}\\right) U_{lith}\n",
"$$\n",
"\n",
"Where $K$ is a fitting parameter, $I$ is the cell current, and $Q_{cell}$ is the cell capacity. To simplify the comparison with the other models, we can rewrite this expression in terms of the hysteresis state variable $h$ given by\n",

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest to use another symbol for $K$ to avoid confusion with rate parameters in other models.

Since the current sigmoid '$K$' is hardcoded in the implementation, quote the value here and perhaps consider the current lower bound (C/1000? C/1e6?) at which the current sigmoid implementation will noticeably no longer obey the intended relation, due to the finite step window.

"\n",
"$$h = 1 - 2 \\, \\text{sigmoid}\\left(-\\frac{KI}{Q_{cell}}\\right).$$\n",
"\n",
"### Hysteresis State Variable (Axen)\n",
"In this model, the state variable $h(z,t)$ is both stoichiometry and time-dependent, and is governed by the following ODEs:\n",
"$$ \\frac{dh(z,t)}{dt} = K_{lith} \\, \\frac{I_{cell}}{Q_{cell}} \\, (1 - h(z,t)) \\qquad \\text{for lithiation}$$\n",

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Even if this is how it's presented in the original literature, it's inherently confusing that the implied range of $h$ is -0.5 to +0.5 here, but -1 to 1 for the general model.

I would recommend that we use this document to collectively paraphrase all the original literature models into a common notation and definition. Do not worry about exactly how the original authors approached symbols or factorisation. The main point of this document should be for PyBaMM users to distinguish the various model options, and understand their basis in scientific literature where appropriate.

"$$ \\frac{dh(z,t)}{dt} = K_{delith} \\, \\frac{I_{cell}}{Q_{cell}} \\, h(z,t) \\qquad \\text{for delithiation}$$\n",
"where $I_{cell}$ is the cell current, $Q_{cell}$ is the cell capacity, and $K_{lith}$ and $K_{delith}$ are the decay rates for lithiation and delithiation, respectively. It is important to highlight that in this form, we need to adjust the sign of the ODE system depending on whether the hysteresis occurs in the anode or in the cathode to align with the sign of the applied current $I_{cell}$.\n",
"\n",
"However, we desire an ODE system with a sign that is invariant to the change in the direction of the applied current. This is achieved by replacing $I_{cell}/Q_{cell}$ with $i_{vol} / (F \\, c_{max} \\, \\epsilon)$ to obtain a general expression that can be applied to both lithiation and delithiation, without the need of adjusting the sign of the formulation$:\n",
"\n",
"$$ \\frac{dh(z,t)}{dt} = -K' \\, \\left(\\frac{1-\\text{sgn}(i_{vol})}{2} + \\text{sgn}(i_{vol}) \\, h(z,t)\\right)$$\n",
"\n",
"$$K' = \\frac{i_{vol}}{F \\, c_{max} \\, \\epsilon} \\, \\left( \\frac{1-\\text{sgn}(i_{vol})}{2} \\, K_{lith} + \\frac{1+\\text{sgn}(i_{vol})}{2} \\, K_{delith}\\right)$$\n",
"\n",
"where $i_{vol}$ is the volumetric interfacial current density of the active material particles showing hysteresis, which is always negative for lithiation and positive for delithiation, regardless of whether these particles are in the anode or in the cathode; $F$ is the Faraday constant, $c_{max}$ is the maximum lithium concentration in the particles and $\\epsilon$ is the active material volume fraction.\n",
"\n",
"The expression to calculate $U$ and the corresponding ODE system is further modified to ensure $h(z,t)$ varies between -1 and 1, thus assuming the same variation range as in the single-state models implemented in PyBaMM.\n",
"\n",
"The resulting ODE governing the hysteresis state variable $h(z,t)$ is:\n",
"\n",
"$$ \\frac{dh(z,t)}{dt} = K' \\, (1 - \\text{sgn}(i_{vol}) \\, h(z,t))$$\n",
"\n",
"### Differential Capacity Hysteresis State (Wycisk)\n",
"As in the Axen model, the state variable $h(z,t)$ is both stoichiometry and time-dependent, and is governed by the following ODE:\n",
"\n",
"$$ \\frac{dh(z,t)}{dt} = \\left(\\frac{k(z) \\cdot I(t)}{Q_{cell}}\\right)\\left(1-\\text{sgn}\\left(\\frac{dz(t)}{dt}\\right) h(z,t)\\right) $$\n",
"\n",
"where $ k(z) $ is given by\n",
"\n",
"$$ k(z) = K \\cdot \\frac{1}{\\left(C_{diff}\\left(z\\right)\\right)^{x}}.$$\n",
"\n",
"Here $C_{diff}(z)$ is the differential capacity with respect to potential, expressed as \n",
"\n",
"$$ C_{diff}(z) = \\frac{dQ}{dU_{eq}(z)},$$\n",
"where $Q$ is the capacity of the phase of active material experiencing the voltage hysteresis and $U_{eq}$ is the average open-circuit potential of the phase of active material experiencing the voltage hysteresis. The remaining parameters are $K$ and $x$ which are both fitting parameters that affect the response of the hysteresis state decay when passing charge in either direction.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparing the different hysteresis state models\n",
"\n",
"We now compare the behavior of the different hysteresis state models. First we define the models, which are all based on the DFN model with a composite negative electrode comprised of two particles phases. One phase is modelled using a single open-circuit potential, while the other is modelled using a hysteresis state model. The positive electrode is modelled using a single open-circuit potential."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"models = {\n",
" \"current sigmoid\": pybamm.lithium_ion.DFN(\n",
" {\n",
" \"open-circuit potential\": ((\"single\", \"current sigmoid\"), \"single\"),\n",
" \"particle phases\": (\"2\", \"1\"),\n",
" \"thermal\": \"lumped\",\n",
" }\n",
" ),\n",
" \"Wycisk\": pybamm.lithium_ion.DFN(\n",
" {\n",
" \"open-circuit potential\": ((\"single\", \"Wycisk\"), \"single\"),\n",
" \"particle phases\": (\"2\", \"1\"),\n",
" \"thermal\": \"lumped\",\n",
" }\n",
" ),\n",
" \"Axen\": pybamm.lithium_ion.DFN(\n",
" {\n",
" \"open-circuit potential\": ((\"single\", \"Axen\"), \"single\"),\n",
" \"particle phases\": (\"2\", \"1\"),\n",
" \"thermal\": \"lumped\",\n",
" }\n",
" ),\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we define the parameter set, which is for a graphite and silicon negative electrode with a NMC positive electrode. We add some "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"parameters = pybamm.ParameterValues(\"Chen2020_composite\")\n",
"parameters.update(\n",
" {\n",
" \"Negative current collector density [kg.m-3]\": 8933.0,\n",
" \"Negative current collector specific heat capacity [J.kg-1.K-1]\": 385.0,\n",
" \"Negative current collector thermal conductivity [W.m-1.K-1]\": 398.0,\n",
" \"Negative electrode density [kg.m-3]\": 1555.0,\n",
" \"Negative electrode specific heat capacity [J.kg-1.K-1]\": 1437.0,\n",
" \"Negative electrode thermal conductivity [W.m-1.K-1]\": 1.58,\n",
" \"Separator density [kg.m-3]\": 1017.0,\n",
" \"Separator specific heat capacity [J.kg-1.K-1]\": 1978.0,\n",
" \"Separator thermal conductivity [W.m-1.K-1]\": 0.34,\n",
" \"Positive electrode density [kg.m-3]\": 3699.0,\n",
" \"Positive electrode specific heat capacity [J.kg-1.K-1]\": 1270.0,\n",
" \"Positive electrode thermal conductivity [W.m-1.K-1]\": 1.04,\n",
" \"Positive current collector density [kg.m-3]\": 2702.0,\n",
" \"Positive current collector specific heat capacity [J.kg-1.K-1]\": 903.0,\n",
" \"Positive current collector thermal conductivity [W.m-1.K-1]\": 238.0,\n",
" \"Total heat transfer coefficient [W.m-2.K-1]\": 5.0,\n",
" },\n",
" check_already_exists=False,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we need to add the additional parameters required by the model. In this example, $K_{lith} = K_{delith} = 100$ for the HS model, but it is recommended to test different values for the decay rate aiming for good agreement to experiment if the information is available, or a smooth transition between the lithiation and delithiation branches of the hysteresis. For example, depending on the case study, an acceptable value for $K_{lith}$ and $K_{delith}$ could be $10$ or lower."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Wycisk parameters\n",
"parameters_wycisk = parameters.copy()\n",
"\n",
"# get the lithiation and delithiation functions and define the average OCP\n",
"lithiation_ocp = parameters_wycisk[\"Secondary: Negative electrode lithiation OCP [V]\"]\n",
"delithiation_ocp = parameters_wycisk[\n",
" \"Secondary: Negative electrode delithiation OCP [V]\"\n",
"]\n",
"\n",
"\n",
"def ocp_avg(sto):\n",
" return (lithiation_ocp(sto) + delithiation_ocp(sto)) / 2\n",
"\n",
"\n",
"parameters_wycisk.update(\n",
" {\n",
" \"Secondary: Negative electrode OCP [V]\": ocp_avg,\n",
" \"Secondary: Negative particle hysteresis decay rate\": 0.005,\n",
" \"Secondary: Negative particle hysteresis switching factor\": 10,\n",
" \"Secondary: Initial hysteresis state in negative electrode\": 0,\n",
" },\n",
" check_already_exists=False,\n",
")\n",
"\n",
"# Axen parameters\n",
"parameters_axen = parameters.copy()\n",
"parameters_axen.update(\n",
" {\n",
" \"Secondary: Negative particle lithiation hysteresis decay rate\": 100,\n",
" \"Secondary: Negative particle delithiation hysteresis decay rate\": 100,\n",
" \"Secondary: Initial hysteresis state in negative electrode\": 0,\n",
" },\n",
" check_already_exists=False,\n",
")\n",
"\n",
"parameter_values = {\n",
" \"current sigmoid\": parameters,\n",
" \"Wycisk\": parameters_wycisk,\n",
" \"Axen\": parameters_axen,\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we define the experiment, which is a discharge and charge cycle at 1C."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"experiment = pybamm.Experiment(\n",
" [\n",
" (\n",
" \"Discharge at 1C for 1 hour or until 2.5 V\",\n",
" \"Rest for 15 minutes\",\n",
" \"Charge at 1C until 4.2 V\",\n",
" \"Hold at 4.2 V until 0.05 C\",\n",
" \"Rest for 15 minutes\",\n",
" ),\n",
" ]\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then loop over the models and solve them."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"solutions = {}\n",
"for name, model in models.items():\n",
" sim = pybamm.Simulation(\n",
" model, experiment=experiment, parameter_values=parameter_values[name]\n",
" )\n",
" solutions[name] = sim.solve(calc_esoh=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we plot the results."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "e65b398609174ea9b8379572db0c36f5",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(FloatSlider(value=0.0, description='t', max=3.003330090763245, step=0.03003330090763245)…"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"<pybamm.plotting.quick_plot.QuickPlot at 0x10603b550>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"output_variables = [\n",
" \"Current [A]\",\n",
" \"Negative electrode secondary stoichiometry\",\n",
" \"Voltage [V]\",\n",
" \"Volume-averaged cell temperature [K]\",\n",
" \"Total heating [W.m-3]\",\n",
" \"Hysteresis electrochemical heating [W.m-3]\",\n",
" \"X-averaged negative electrode secondary hysteresis state\",\n",
" \"X-averaged negative electrode secondary open-circuit potential [V]\",\n",
"]\n",
"\n",
"pybamm.dynamic_plot(\n",
" list(solutions.values()),\n",
" labels=list(solutions.keys()),\n",
" colors=[\"black\", \"blue\", \"red\"],\n",
" linestyles=[\"-\", \"--\", \"-.\"],\n",
" output_variables=output_variables,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
from .base_ocp import BaseOpenCircuitPotential
from .base_hysteresis_ocp import BaseHysteresisOpenCircuitPotential
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', 'axen_ocp']
__all__ = ['base_ocp', 'base_hysteresis_ocp', 'current_sigmoid_ocp', 'msmr_ocp', 'single_ocp', 'wycisk_ocp', 'axen_ocp']
Loading
Loading