diff --git a/CHANGELOG.md b/CHANGELOG.md index 7df41530..b3ce61ee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#452](https://github.com/pybop-team/PyBOP/issues/452) - Extends `cell_mass` and `approximate_capacity` for half-cell models. - [#544](https://github.com/pybop-team/PyBOP/issues/544) - Allows iterative plotting using `StandardPlot`. - [#541](https://github.com/pybop-team/PyBOP/pull/541) - Adds `ScaledLogLikelihood` and `BaseMetaLikelihood` classes. - [#409](https://github.com/pybop-team/PyBOP/pull/409) - Adds plotting and convergence methods for Monte Carlo sampling. Includes open-access Tesla 4680 dataset for Bayesian inference example. Fixes transformations for sampling. diff --git a/examples/notebooks/energy_based_electrode_design.ipynb b/examples/notebooks/energy_based_electrode_design.ipynb index 50f8ed5f..96871b9f 100644 --- a/examples/notebooks/energy_based_electrode_design.ipynb +++ b/examples/notebooks/energy_based_electrode_design.ipynb @@ -76,6 +76,7 @@ "outputs": [], "source": [ "import numpy as np\n", + "from pybamm import Parameter\n", "\n", "import pybop\n", "\n", @@ -125,6 +126,24 @@ "outputs": [], "source": [ "parameter_set = pybop.ParameterSet.pybamm(\"Chen2020\")\n", + "parameter_set.update(\n", + " {\n", + " \"Electrolyte density [kg.m-3]\": Parameter(\"Separator density [kg.m-3]\"),\n", + " \"Negative electrode active material density [kg.m-3]\": Parameter(\n", + " \"Negative electrode density [kg.m-3]\"\n", + " ),\n", + " \"Negative electrode carbon-binder density [kg.m-3]\": Parameter(\n", + " \"Negative electrode density [kg.m-3]\"\n", + " ),\n", + " \"Positive electrode active material density [kg.m-3]\": Parameter(\n", + " \"Positive electrode density [kg.m-3]\"\n", + " ),\n", + " \"Positive electrode carbon-binder density [kg.m-3]\": Parameter(\n", + " \"Positive electrode density [kg.m-3]\"\n", + " ),\n", + " },\n", + " check_already_exists=False,\n", + ")\n", "model = pybop.lithium_ion.SPMe(parameter_set=parameter_set)" ] }, diff --git a/examples/scripts/linked_parameters.py b/examples/scripts/linked_parameters.py index d368d001..9d4797bd 100644 --- a/examples/scripts/linked_parameters.py +++ b/examples/scripts/linked_parameters.py @@ -1,18 +1,40 @@ -import pybamm +from pybamm import Parameter import pybop # The aim of this script is to show how to systematically update # design parameters which depend on the optimisation parameters. -# Define parameter set and model +# Define parameter set and additional parameters needed for the cost function parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True) +parameter_set.update( + { + "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"), + "Negative electrode active material density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Negative electrode carbon-binder density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Positive electrode active material density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode carbon-binder density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + }, + check_already_exists=False, +) + +# Link parameters parameter_set.update( { "Positive electrode porosity": 1 - - pybamm.Parameter("Positive electrode active material volume fraction") + - Parameter("Positive electrode active material volume fraction") } ) + +# Define model model = pybop.lithium_ion.SPMe(parameter_set=parameter_set) # Fitting parameters diff --git a/examples/scripts/maximising_energy.py b/examples/scripts/maximising_energy.py index 6d0a82d9..23364664 100644 --- a/examples/scripts/maximising_energy.py +++ b/examples/scripts/maximising_energy.py @@ -1,3 +1,5 @@ +from pybamm import Parameter + import pybop # A design optimisation example loosely based on work by L.D. Couto @@ -9,8 +11,28 @@ # electrode widths, particle radii, volume fractions and # separator width. -# Define parameter set and model +# Define parameter set and additional parameters needed for the cost function parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True) +parameter_set.update( + { + "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"), + "Negative electrode active material density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Negative electrode carbon-binder density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Positive electrode active material density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode carbon-binder density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + }, + check_already_exists=False, +) + +# Define model model = pybop.lithium_ion.SPMe(parameter_set=parameter_set) # Fitting parameters diff --git a/examples/scripts/maximising_power.py b/examples/scripts/maximising_power.py index c72a07f1..d1a8f857 100644 --- a/examples/scripts/maximising_power.py +++ b/examples/scripts/maximising_power.py @@ -1,7 +1,29 @@ +from pybamm import Parameter + import pybop -# Define parameter set and model +# Define parameter set and additional parameters needed for the cost function parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True) +parameter_set.update( + { + "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"), + "Negative electrode active material density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Negative electrode carbon-binder density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Positive electrode active material density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode carbon-binder density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + }, + check_already_exists=False, +) + +# Define model model = pybop.lithium_ion.SPMe(parameter_set=parameter_set) # Define useful quantities diff --git a/pybop/models/lithium_ion/base_echem.py b/pybop/models/lithium_ion/base_echem.py index 683188fb..ade7e8fa 100644 --- a/pybop/models/lithium_ion/base_echem.py +++ b/pybop/models/lithium_ion/base_echem.py @@ -2,6 +2,7 @@ import warnings from typing import Optional +from pybamm import LithiumIonParameters from pybamm import lithium_ion as pybamm_lithium_ion from pybop.models.base_model import BaseModel, Inputs @@ -86,6 +87,7 @@ def __init__( self._disc = None self._electrode_soh = pybamm_lithium_ion.electrode_soh + self._electrode_soh_half_cell = pybamm_lithium_ion.electrode_soh_half_cell self.geometric_parameters = self.set_geometric_parameters() def _check_params( @@ -213,30 +215,36 @@ def cell_mass(self, parameter_set: Optional[ParameterSet] = None): parameter_set = parameter_set or self._parameter_set def mass_density( - active_material_vol_frac, density, porosity, electrolyte_density + active_material_vol_frac, + density, + porosity, + electrolyte_density, + carbon_binder_domain_density, ): - return (active_material_vol_frac * density) + ( - porosity * electrolyte_density + return ( + (active_material_vol_frac * density) + + (porosity * electrolyte_density) + + (1.0 - active_material_vol_frac - porosity) + * carbon_binder_domain_density ) def area_density(thickness, mass_density): return thickness * mass_density - # Approximations due to SPM(e) parameter set limitations - electrolyte_density = parameter_set["Separator density [kg.m-3]"] - # Calculate mass densities positive_mass_density = mass_density( parameter_set["Positive electrode active material volume fraction"], - parameter_set["Positive electrode density [kg.m-3]"], + parameter_set["Positive electrode active material density [kg.m-3]"], parameter_set["Positive electrode porosity"], - electrolyte_density, + parameter_set["Electrolyte density [kg.m-3]"], + parameter_set["Positive electrode carbon-binder density [kg.m-3]"], ) negative_mass_density = mass_density( parameter_set["Negative electrode active material volume fraction"], - parameter_set["Negative electrode density [kg.m-3]"], + parameter_set["Negative electrode active material density [kg.m-3]"], parameter_set["Negative electrode porosity"], - electrolyte_density, + parameter_set["Electrolyte density [kg.m-3]"], + parameter_set["Negative electrode carbon-binder density [kg.m-3]"], ) # Calculate area densities @@ -248,7 +256,7 @@ def area_density(thickness, mass_density): ) separator_area_density = area_density( parameter_set["Separator thickness [m]"], - parameter_set["Separator porosity"] * electrolyte_density, + parameter_set["Separator density [kg.m-3]"], ) positive_cc_area_density = area_density( parameter_set["Positive current collector thickness [m]"], @@ -279,8 +287,8 @@ def area_density(thickness, mass_density): def approximate_capacity(self, parameter_set: Optional[ParameterSet] = None): """ Calculate an estimate for the nominal cell capacity. The estimate is computed - by dividing the theoretical energy (in watt-hours) by the average open circuit - potential (voltage) of the cell. + by estimating the capacity of the positive electrode that lies between the + stoichiometric limits corresponding to the upper and lower voltage limits. Parameters ---------- @@ -290,39 +298,31 @@ def approximate_capacity(self, parameter_set: Optional[ParameterSet] = None): Returns ------- float - The estimate of the nominal cell capacity. + The estimate of the nominal cell capacity [A.h]. """ parameter_set = parameter_set or self._parameter_set - # Calculate theoretical energy density - theoretical_energy = self._electrode_soh.calculate_theoretical_energy( - parameter_set - ) - - # Extract stoichiometries and compute mean values - ( - min_sto_neg, - max_sto_neg, - min_sto_pos, - max_sto_pos, - ) = self._electrode_soh.get_min_max_stoichiometries(parameter_set) - mean_sto_neg = (min_sto_neg + max_sto_neg) / 2 - mean_sto_pos = (min_sto_pos + max_sto_pos) / 2 - - # Calculate average voltage - positive_electrode_ocp = parameter_set["Positive electrode OCP [V]"] - negative_electrode_ocp = parameter_set["Negative electrode OCP [V]"] - try: - average_voltage = positive_electrode_ocp( - mean_sto_pos - ) - negative_electrode_ocp(mean_sto_neg) - except Exception as e: - raise ValueError(f"Error in average voltage calculation: {e}") from e - - # Calculate the capacity estimate - approximate_capacity = theoretical_energy / average_voltage - - return ParameterSet.evaluate_symbol(approximate_capacity, parameter_set) + # Calculate the theoretical capacity in the limit of low current + if self.pybamm_model.options["working electrode"] == "positive": + ( + max_sto_p, + min_sto_p, + ) = self._electrode_soh_half_cell.get_min_max_stoichiometries(parameter_set) + else: + ( + min_sto_n, + max_sto_n, + min_sto_p, + max_sto_p, + ) = self._electrode_soh.get_min_max_stoichiometries(parameter_set) + # Note that the stoichiometric limits correspond to 0 and 100% SOC. + # Stoichiometric balancing is performed within get_min_max_stoichiometries + # such that the capacity accessible between the limits should be the same + # for both electrodes, so we consider just the positive electrode below. + + Q_p = LithiumIonParameters().p.prim.Q_init + theoretical_capacity = Q_p * (max_sto_p - min_sto_p) + return ParameterSet.evaluate_symbol(theoretical_capacity, parameter_set) def set_geometric_parameters(self): """ diff --git a/tests/integration/test_half_cell_model.py b/tests/integration/test_half_cell_model.py new file mode 100644 index 00000000..f3c3a270 --- /dev/null +++ b/tests/integration/test_half_cell_model.py @@ -0,0 +1,163 @@ +import numpy as np +import pytest +from pybamm import Parameter + +import pybop + + +class TestHalfCellModel: + """ + A class to test optimisation of a PyBaMM half-cell model. + """ + + @pytest.fixture(autouse=True) + def setup(self): + self.sigma0 = 0.002 + self.ground_truth = np.clip( + np.asarray([0.5]) + np.random.normal(loc=0.0, scale=0.05, size=1), + a_min=0.4, + a_max=0.75, + ) + + @pytest.fixture + def model(self): + options = {"working electrode": "positive"} + parameter_set = pybop.lithium_ion.SPM(options=options).default_parameter_values + parameter_set.update( + { + "Electrolyte density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Negative current collector density [kg.m-3]": 0.0, + "Negative current collector thickness [m]": 0.0, + "Negative electrode active material density [kg.m-3]": 0.0, + "Negative electrode active material volume fraction": 0.0, + "Negative electrode porosity": 0.0, + "Negative electrode carbon-binder density [kg.m-3]": 0.0, + "Positive current collector density [kg.m-3]": 0.0, + "Positive current collector thickness [m]": 0.0, + "Positive electrode active material density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode carbon-binder density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode density [kg.m-3]": 3262.0, + "Separator density [kg.m-3]": 0.0, + }, + check_already_exists=False, + ) + x = self.ground_truth + parameter_set.update( + { + "Positive electrode active material volume fraction": x[0], + } + ) + return pybop.lithium_ion.SPM(parameter_set=parameter_set, options=options) + + @pytest.fixture + def parameters(self): + return pybop.Parameters( + pybop.Parameter( + "Positive electrode active material volume fraction", + prior=pybop.Uniform(0.4, 0.75), + # no bounds + ), + ) + + @pytest.fixture(params=[0.4]) + def init_soc(self, request): + return request.param + + def noise(self, sigma, values): + return np.random.normal(0, sigma, values) + + @pytest.fixture + def fitting_cost(self, model, parameters, init_soc): + # Form dataset + solution = self.get_data(model, init_soc) + dataset = pybop.Dataset( + { + "Time [s]": solution["Time [s]"].data, + "Current function [A]": solution["Current [A]"].data, + "Voltage [V]": solution["Voltage [V]"].data + + self.noise(self.sigma0, len(solution["Time [s]"].data)), + } + ) + + # Define the cost to optimise + problem = pybop.FittingProblem(model, parameters, dataset) + return pybop.SumSquaredError(problem) + + @pytest.mark.integration + def test_fitting_costs(self, fitting_cost): + x0 = fitting_cost.parameters.initial_value() + optim = pybop.CuckooSearch( + cost=fitting_cost, + sigma0=0.03, + max_iterations=250, + max_unchanged_iterations=35, + ) + + initial_cost = optim.cost(optim.parameters.initial_value()) + results = optim.run() + + # Assertions + if not np.allclose(x0, self.ground_truth, atol=1e-5): + if optim.minimising: + assert initial_cost > results.final_cost + else: + assert initial_cost < results.final_cost + np.testing.assert_allclose(results.x, self.ground_truth, atol=1.5e-2) + + @pytest.fixture + def design_cost(self, model): + initial_state = {"Initial SoC": 1.0} + parameters = pybop.Parameters( + pybop.Parameter( + "Positive electrode thickness [m]", + prior=pybop.Gaussian(5e-05, 5e-06), + bounds=[2e-06, 10e-05], + ), + ) + experiment = pybop.Experiment( + ["Discharge at 1C until 3.5 V (5 seconds period)"], + ) + + problem = pybop.DesignProblem( + model, + parameters, + experiment=experiment, + initial_state=initial_state, + update_capacity=True, + ) + return pybop.GravimetricEnergyDensity(problem) + + @pytest.mark.integration + def test_design_costs(self, design_cost): + optim = pybop.CuckooSearch( + design_cost, + max_iterations=15, + allow_infeasible_solutions=False, + ) + initial_values = optim.parameters.initial_value() + initial_cost = optim.cost(initial_values) + results = optim.run() + + # Assertions + assert initial_cost < results.final_cost + for i, _ in enumerate(results.x): + assert results.x[i] < initial_values[i] + + def get_data(self, model, init_soc): + initial_state = {"Initial SoC": init_soc} + experiment = pybop.Experiment( + [ + ( + "Discharge at 0.5C for 3 minutes (4 second period)", + "Charge at 0.5C for 3 minutes (4 second period)", + ), + ] + ) + sim = model.predict(initial_state=initial_state, experiment=experiment) + return sim diff --git a/tests/integration/test_weighted_cost.py b/tests/integration/test_weighted_cost.py index bd5baca9..13df11cd 100644 --- a/tests/integration/test_weighted_cost.py +++ b/tests/integration/test_weighted_cost.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from pybamm import Parameter import pybop @@ -21,6 +22,24 @@ def setup(self): @pytest.fixture def model(self): parameter_set = pybop.ParameterSet.pybamm("Chen2020") + parameter_set.update( + { + "Electrolyte density [kg.m-3]": Parameter("Separator density [kg.m-3]"), + "Negative electrode active material density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Negative electrode carbon-binder density [kg.m-3]": Parameter( + "Negative electrode density [kg.m-3]" + ), + "Positive electrode active material density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode carbon-binder density [kg.m-3]": Parameter( + "Positive electrode density [kg.m-3]" + ), + }, + check_already_exists=False, + ) x = self.ground_truth parameter_set.update( { diff --git a/tests/unit/test_cost.py b/tests/unit/test_cost.py index a8247c5a..e9ad01dc 100644 --- a/tests/unit/test_cost.py +++ b/tests/unit/test_cost.py @@ -204,7 +204,28 @@ def test_sumofpower(self, problem): @pytest.fixture def design_problem(self, parameters, experiment, signal): - model = pybop.lithium_ion.SPM() + parameter_set = pybop.ParameterSet.pybamm("Chen2020") + parameter_set.update( + { + "Electrolyte density [kg.m-3]": pybamm.Parameter( + "Separator density [kg.m-3]" + ), + "Negative electrode active material density [kg.m-3]": pybamm.Parameter( + "Negative electrode density [kg.m-3]" + ), + "Negative electrode carbon-binder density [kg.m-3]": pybamm.Parameter( + "Negative electrode density [kg.m-3]" + ), + "Positive electrode active material density [kg.m-3]": pybamm.Parameter( + "Positive electrode density [kg.m-3]" + ), + "Positive electrode carbon-binder density [kg.m-3]": pybamm.Parameter( + "Positive electrode density [kg.m-3]" + ), + }, + check_already_exists=False, + ) + model = pybop.lithium_ion.SPM(parameter_set=parameter_set) return pybop.DesignProblem( model, parameters,