diff --git a/CHANGELOG.md b/CHANGELOG.md index e879450e..6f4ed5b7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#249](https://github.com/pybop-team/PyBOP/pull/249) - Add WeppnerHuggins model and GITT example. - [#304](https://github.com/pybop-team/PyBOP/pull/304) - Decreases the testing suite completion time. - [#301](https://github.com/pybop-team/PyBOP/pull/301) - Updates default echem solver to "fast with events" mode. - [#251](https://github.com/pybop-team/PyBOP/pull/251) - Increment PyBaMM > v23.5, remove redundant tests within integration tests, increment citation version, fix examples with incorrect model definitions. diff --git a/examples/scripts/gitt.py b/examples/scripts/gitt.py new file mode 100644 index 00000000..d31ac2d4 --- /dev/null +++ b/examples/scripts/gitt.py @@ -0,0 +1,70 @@ +import numpy as np + +import pybop + +# Define model +parameter_set = pybop.ParameterSet.pybamm("Xu2019") +model = pybop.lithium_ion.SPM( + parameter_set=parameter_set, options={"working electrode": "positive"} +) + +# Generate data +sigma = 0.005 +t_eval = np.arange(0, 150, 2) +values = model.predict(t_eval=t_eval) +corrupt_values = values["Voltage [V]"].data + np.random.normal(0, sigma, len(t_eval)) + +# Form dataset +dataset = pybop.Dataset( + { + "Time [s]": t_eval, + "Current function [A]": values["Current [A]"].data, + "Voltage [V]": corrupt_values, + } +) + +# Define parameter set +parameter_set.update( + { + "Reference OCP [V]": 4.1821, + "Derivative of the OCP wrt stoichiometry [V]": -1.38636, + }, + check_already_exists=False, +) + +# Define the cost to optimise +model = pybop.lithium_ion.WeppnerHuggins(parameter_set=parameter_set) + +parameters = [ + pybop.Parameter( + "Positive electrode diffusivity [m2.s-1]", + prior=pybop.Gaussian(5e-14, 1e-13), + bounds=[1e-16, 1e-11], + true_value=parameter_set["Positive electrode diffusivity [m2.s-1]"], + ), +] + +problem = pybop.FittingProblem( + model, + parameters, + dataset, + signal=["Voltage [V]"], +) + +cost = pybop.RootMeanSquaredError(problem) + +# Build the optimisation problem +optim = pybop.Optimisation(cost=cost, optimiser=pybop.PSO, verbose=True) + +# Run the optimisation problem +x, final_cost = optim.run() +print("Estimated parameters:", x) + +# Plot the timeseries output +pybop.quick_plot(problem, parameter_values=x, title="Optimised Comparison") + +# Plot convergence +pybop.plot_convergence(optim) + +# Plot the parameter traces +pybop.plot_parameters(optim) diff --git a/pybop/__init__.py b/pybop/__init__.py index c25abdc7..034dcb4a 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -46,8 +46,7 @@ from ._utils import is_numeric # -# Cost class -# Problem class +# Problem classes # from .problems.base_problem import BaseProblem from .problems.fitting_problem import FittingProblem diff --git a/pybop/models/lithium_ion/__init__.py b/pybop/models/lithium_ion/__init__.py index 80e2d2cc..6e85f674 100644 --- a/pybop/models/lithium_ion/__init__.py +++ b/pybop/models/lithium_ion/__init__.py @@ -2,4 +2,4 @@ # Import lithium ion based models # from .base_echem import EChemBaseModel -from .echem import SPM, SPMe, DFN, MPM, MSMR +from .echem import SPM, SPMe, DFN, MPM, MSMR, WeppnerHuggins diff --git a/pybop/models/lithium_ion/base_echem.py b/pybop/models/lithium_ion/base_echem.py index f681bd89..d47a6c6f 100644 --- a/pybop/models/lithium_ion/base_echem.py +++ b/pybop/models/lithium_ion/base_echem.py @@ -73,16 +73,24 @@ def _check_params( """ parameter_set = parameter_set or self._parameter_set - electrode_params = [ - ( - "Negative electrode active material volume fraction", - "Negative electrode porosity", - ), - ( - "Positive electrode active material volume fraction", - "Positive electrode porosity", - ), - ] + if self.pybamm_model.options["working electrode"] == "positive": + electrode_params = [ + ( + "Positive electrode active material volume fraction", + "Positive electrode porosity", + ), + ] + else: + electrode_params = [ + ( + "Negative electrode active material volume fraction", + "Negative electrode porosity", + ), + ( + "Positive electrode active material volume fraction", + "Positive electrode porosity", + ), + ] related_parameters = { key: inputs.get(key) if inputs and key in inputs else parameter_set[key] diff --git a/pybop/models/lithium_ion/echem.py b/pybop/models/lithium_ion/echem.py index 30e9e6d0..ba4d0140 100644 --- a/pybop/models/lithium_ion/echem.py +++ b/pybop/models/lithium_ion/echem.py @@ -1,6 +1,7 @@ import pybamm from .base_echem import EChemBaseModel +from .weppner_huggins import BaseWeppnerHuggins class SPM(EChemBaseModel): @@ -270,3 +271,50 @@ def __init__( spatial_methods=spatial_methods, solver=solver, ) + + +class WeppnerHuggins(EChemBaseModel): + """ + Represents the Weppner & Huggins model to fit diffusion coefficients to GITT data. + + Parameters + ---------- + name: str, optional + A name for the model instance, defaults to "Weppner & Huggins model". + parameter_set: pybamm.ParameterValues or dict, optional + A dictionary or a ParameterValues object containing the parameters for the model. If None, the default parameters are used. + geometry: dict, optional + A dictionary defining the model's geometry. If None, the default geometry is used. + submesh_types: dict, optional + A dictionary defining the types of submeshes to use. If None, the default submesh types are used. + var_pts: dict, optional + A dictionary specifying the number of points for each variable for discretization. If None, the default variable points are used. + spatial_methods: dict, optional + A dictionary specifying the spatial methods for discretization. If None, the default spatial methods are used. + solver: pybamm.Solver, optional + The solver to use for simulating the model. If None, the default solver is used. + """ + + def __init__( + self, + name="Weppner & Huggins model", + parameter_set=None, + geometry=None, + submesh_types=None, + var_pts=None, + spatial_methods=None, + solver=None, + ): + self.pybamm_model = BaseWeppnerHuggins() + self._unprocessed_model = self.pybamm_model + + super().__init__( + model=self.pybamm_model, + name=name, + parameter_set=parameter_set, + geometry=geometry, + submesh_types=submesh_types, + var_pts=var_pts, + spatial_methods=spatial_methods, + solver=solver, + ) diff --git a/pybop/models/lithium_ion/weppner_huggins.py b/pybop/models/lithium_ion/weppner_huggins.py new file mode 100644 index 00000000..d656ea72 --- /dev/null +++ b/pybop/models/lithium_ion/weppner_huggins.py @@ -0,0 +1,113 @@ +# +# Weppner Huggins Model +# +import numpy as np +import pybamm + + +class BaseWeppnerHuggins(pybamm.lithium_ion.BaseModel): + """WeppnerHuggins Model for GITT. Credit: pybamm-param team. + + Parameters + ---------- + name : str, optional + The name of the model. + """ + + def __init__(self, name="Weppner & Huggins model"): + super().__init__({}, name) + + pybamm.citations.register(""" + @article{Weppner1977, + title={{Determination of the kinetic parameters + of mixed-conducting electrodes and application to the system Li3Sb}}, + author={Weppner, W and Huggins, R A}, + journal={Journal of The Electrochemical Society}, + volume={124}, + number={10}, + pages={1569}, + year={1977}, + publisher={IOP Publishing} + } + """) + + # `self.param` is a class containing all the relevant parameters and functions for + # this model. These are purely symbolic at this stage, and will be set by the + # `ParameterValues` class when the model is processed. + self.options["working electrode"] = "positive" + self._summary_variables = [] + + t = pybamm.t + ###################### + # Parameters + ###################### + + d_s = pybamm.Parameter("Positive electrode diffusivity [m2.s-1]") + + c_s_max = pybamm.Parameter( + "Maximum concentration in positive electrode [mol.m-3]" + ) + + i_app = self.param.current_density_with_time + + U = pybamm.Parameter("Reference OCP [V]") + + U_prime = pybamm.Parameter("Derivative of the OCP wrt stoichiometry [V]") + + epsilon = pybamm.Parameter("Positive electrode active material volume fraction") + + r_particle = pybamm.Parameter("Positive particle radius [m]") + + a = 3 * (epsilon / r_particle) + + l_w = self.param.p.L + + ###################### + # Governing equations + ###################### + u_surf = ( + (2 / (np.pi**0.5)) + * (i_app / ((d_s**0.5) * a * self.param.F * l_w)) + * (t**0.5) + ) + # Linearised voltage + V = U + (U_prime * u_surf) / c_s_max + ###################### + # (Some) variables + ###################### + self.variables = { + "Voltage [V]": V, + "Time [s]": t, + } + + @property + def default_geometry(self): + return {} + + @property + def default_parameter_values(self): + parameter_values = pybamm.ParameterValues("Xu2019") + parameter_values.update( + { + "Reference OCP [V]": 4.1821, + "Derivative of the OCP wrt stoichiometry [V]": -1.38636, + }, + check_already_exists=False, + ) + return parameter_values + + @property + def default_submesh_types(self): + return {} + + @property + def default_var_pts(self): + return {} + + @property + def default_spatial_methods(self): + return {} + + @property + def default_solver(self): + return pybamm.DummySolver() diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index ba545a81..3137f6f2 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -18,6 +18,7 @@ class TestModels: pybop.lithium_ion.DFN(), pybop.lithium_ion.MPM(), pybop.lithium_ion.MSMR(options={"number of MSMR reactions": ("6", "4")}), + pybop.lithium_ion.WeppnerHuggins(), pybop.empirical.Thevenin(), ] ) @@ -61,10 +62,15 @@ def test_predict_with_inputs(self, model): # Define inputs t_eval = np.linspace(0, 10, 100) if isinstance(model, (pybop.lithium_ion.EChemBaseModel)): - inputs = { - "Negative electrode active material volume fraction": 0.52, - "Positive electrode active material volume fraction": 0.63, - } + if model.pybamm_model.options["working electrode"] == "positive": + inputs = { + "Positive electrode active material volume fraction": 0.63, + } + else: + inputs = { + "Negative electrode active material volume fraction": 0.52, + "Positive electrode active material volume fraction": 0.63, + } elif isinstance(model, (pybop.empirical.Thevenin)): inputs = { "R0 [Ohm]": 0.0002,