diff --git a/CHANGELOG.md b/CHANGELOG.md index 577168ed..e038233d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#532](https://github.com/pybop-team/PyBOP/issues/532) - Adds `linked_parameters` example script which shows how to update linked parameters during design optimisation. - [#529](https://github.com/pybop-team/PyBOP/issues/529) - Adds `GravimetricPowerDensity` and `VolumetricPowerDensity` costs, along with the mathjax extension for Sphinx. ## Optimisations diff --git a/README.md b/README.md index 9bfcc1e4..418405b7 100644 --- a/README.md +++ b/README.md @@ -64,7 +64,7 @@ These include a wide variety of optimisation problems that require careful consi Explore our [example notebooks](https://github.com/pybop-team/PyBOP/blob/develop/examples) for hands-on demonstrations: -- [Gravimetric design optimisation (SPM)](https://github.com/pybop-team/PyBOP/blob/develop/examples/notebooks/spm_electrode_design.ipynb) +- [Gravimetric design optimisation (SPM)](https://github.com/pybop-team/PyBOP/blob/develop/examples/notebooks/energy_based_electrode_design.ipynb) - [Non-linear constrained ECM parameter identification](https://github.com/pybop-team/PyBOP/blob/develop/examples/notebooks/ecm_trust-constr.ipynb) - [Optimiser comparison for parameter identification](https://github.com/pybop-team/PyBOP/blob/develop/examples/notebooks/multi_optimiser_identification.ipynb) - [Parameter identification for spatial pouch cell model](https://github.com/pybop-team/PyBOP/blob/develop/examples/notebooks/pouch_cell_identification.ipynb) diff --git a/examples/scripts/linked_parameters.py b/examples/scripts/linked_parameters.py new file mode 100644 index 00000000..3e25490d --- /dev/null +++ b/examples/scripts/linked_parameters.py @@ -0,0 +1,67 @@ +import pybamm + +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 +parameter_set = pybop.ParameterSet.pybamm("Chen2020", formation_concentrations=True) +parameter_set.update( + { + "Positive electrode porosity": 1 + - pybamm.Parameter("Positive electrode active material volume fraction") + } +) +model = pybop.lithium_ion.SPMe(parameter_set=parameter_set) + +# Fitting parameters +parameters = pybop.Parameters( + pybop.Parameter( + "Positive electrode thickness [m]", + prior=pybop.Gaussian(7.56e-05, 0.1e-05), + bounds=[65e-06, 10e-05], + ), + pybop.Parameter( + "Positive electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.15), + bounds=[0.1, 0.9], + ), +) + +# Define test protocol +experiment = pybop.Experiment( + [ + "Discharge at 1C until 2.5 V (5 seconds period)", + "Hold at 2.5 V for 30 minutes or until 10 mA (5 seconds period)", + ], +) +signal = ["Voltage [V]", "Current [A]"] + +# Generate problem +problem = pybop.DesignProblem( + model, + parameters, + experiment, + signal=signal, + initial_state={"Initial SoC": 1.0}, + update_capacity=True, +) + +# Define the cost +cost = pybop.GravimetricEnergyDensity(problem) + +# Run optimisation +optim = pybop.XNES( + cost, verbose=True, allow_infeasible_solutions=False, max_iterations=10 +) +results = optim.run() +print("Estimated parameters:", results.x) +print(f"Initial gravimetric energy density: {cost(optim.x0):.2f} Wh.kg-1") +print(f"Optimised gravimetric energy density: {cost(results.x):.2f} Wh.kg-1") + +# Plot the timeseries output +pybop.quick_plot(problem, problem_inputs=results.x, title="Optimised Comparison") + +# Plot the cost landscape with optimisation path +pybop.plot2d(optim, steps=5) diff --git a/pybop/models/lithium_ion/base_echem.py b/pybop/models/lithium_ion/base_echem.py index d42b8078..4e7f3a7f 100644 --- a/pybop/models/lithium_ion/base_echem.py +++ b/pybop/models/lithium_ion/base_echem.py @@ -1,3 +1,4 @@ +import sys import warnings from typing import Optional @@ -136,9 +137,12 @@ def _check_params( } for material_vol_fraction, porosity in electrode_params: - if ( + total_vol_fraction = ( related_parameters[material_vol_fraction] + related_parameters[porosity] - > 1 + ) + if ( + ParameterSet.evaluate_symbol(total_vol_fraction, parameter_set) + > 1 + sys.float_info.epsilon ): if self.param_check_counter <= len(electrode_params): infeasibility_warning = "Non-physical point encountered - [{material_vol_fraction} + {porosity}] > 1.0!" @@ -182,8 +186,10 @@ def cell_volume(self, parameter_set: Optional[ParameterSet] = None): parameter_set["Electrode height [m]"] * parameter_set["Electrode width [m]"] ) - # Calculate and return total cell volume - return cross_sectional_area * cell_thickness + # Calculate total cell volume + cell_volume = cross_sectional_area * cell_thickness + + return ParameterSet.evaluate_symbol(cell_volume, parameter_set) def cell_mass(self, parameter_set: Optional[ParameterSet] = None): """ @@ -258,7 +264,7 @@ def area_density(thickness, mass_density): parameter_set["Electrode height [m]"] * parameter_set["Electrode width [m]"] ) - # Calculate and return total cell mass + # Calculate total cell mass total_area_density = ( positive_area_density + negative_area_density @@ -266,7 +272,9 @@ def area_density(thickness, mass_density): + positive_cc_area_density + negative_cc_area_density ) - return cross_sectional_area * total_area_density + cell_mass = cross_sectional_area * total_area_density + + return ParameterSet.evaluate_symbol(cell_mass, parameter_set) def approximate_capacity(self, parameter_set: Optional[ParameterSet] = None): """ @@ -314,7 +322,7 @@ def approximate_capacity(self, parameter_set: Optional[ParameterSet] = None): # Calculate the capacity estimate approximate_capacity = theoretical_energy / average_voltage - return approximate_capacity + return ParameterSet.evaluate_symbol(approximate_capacity, parameter_set) def set_geometric_parameters(self): """ diff --git a/pybop/parameters/parameter_set.py b/pybop/parameters/parameter_set.py index 8cfc2c02..4c4e5347 100644 --- a/pybop/parameters/parameter_set.py +++ b/pybop/parameters/parameter_set.py @@ -1,7 +1,18 @@ import json import types +from numbers import Number +from typing import Union -from pybamm import LithiumIonParameters, ParameterValues, parameter_sets +import numpy as np +from pybamm import ( + FunctionParameter, + LithiumIonParameters, + Parameter, + ParameterValues, + Scalar, + Symbol, + parameter_sets, +) class ParameterSet: @@ -34,6 +45,33 @@ def __setitem__(self, key, value): def __getitem__(self, key): return self.params[key] + @staticmethod + def evaluate_symbol(symbol: Union[Symbol, Number], params: dict): + """ + Evaluate a parameter in the parameter set. + + Parameters + ---------- + symbol : pybamm.Symbol or Number + The parameter to evaluate. + + Returns + ------- + float + The value of the parameter. + """ + if isinstance(symbol, (Number, np.float64)): + return symbol + if isinstance(symbol, Scalar): + return symbol.value + if isinstance(symbol, (Parameter, FunctionParameter)): + return ParameterSet.evaluate_symbol(params[symbol.name], params) + new_children = [ + Scalar(ParameterSet.evaluate_symbol(child, params)) + for child in symbol.children + ] + return symbol.create_copy(new_children).evaluate() + def keys(self) -> list: """ A list of parameter names diff --git a/pybop/problems/fitting_problem.py b/pybop/problems/fitting_problem.py index 25609b6b..16a1f478 100644 --- a/pybop/problems/fitting_problem.py +++ b/pybop/problems/fitting_problem.py @@ -101,11 +101,7 @@ def set_initial_state(self, initial_state: Optional[dict] = None): self.initial_state = initial_state - def evaluate( - self, - inputs: Inputs, - update_capacity=False, - ) -> dict[str, np.ndarray[np.float64]]: + def evaluate(self, inputs: Inputs) -> dict[str, np.ndarray[np.float64]]: """ Evaluate the model with the given parameters and return the signal. @@ -121,7 +117,7 @@ def evaluate( """ inputs = self.parameters.verify(inputs) if self.eis: - return self._evaluateEIS(inputs, update_capacity=update_capacity) + return self._evaluateEIS(inputs) else: try: sol = self._model.simulate( @@ -139,9 +135,7 @@ def evaluate( for signal in (self.signal + self.additional_variables) } - def _evaluateEIS( - self, inputs: Inputs, update_capacity=False - ) -> dict[str, np.ndarray[np.float64]]: + def _evaluateEIS(self, inputs: Inputs) -> dict[str, np.ndarray[np.float64]]: """ Evaluate the model with the given parameters and return the signal. diff --git a/pybop/problems/multi_fitting_problem.py b/pybop/problems/multi_fitting_problem.py index d0bb997e..72dd70c3 100644 --- a/pybop/problems/multi_fitting_problem.py +++ b/pybop/problems/multi_fitting_problem.py @@ -76,7 +76,7 @@ def set_initial_state(self, initial_state: Optional[dict] = None): for problem in self.problems: problem.set_initial_state(initial_state) - def evaluate(self, inputs: Inputs, eis=False): + def evaluate(self, inputs: Inputs): """ Evaluate the model with the given parameters and return the signal. diff --git a/tests/unit/test_parameter_sets.py b/tests/unit/test_parameter_sets.py index 347dfde9..3c4c116b 100644 --- a/tests/unit/test_parameter_sets.py +++ b/tests/unit/test_parameter_sets.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from pybamm import FunctionParameter, Parameter, Scalar import pybop @@ -136,3 +137,19 @@ def test_set_formation_concentrations(self): assert ( parameter_set["Initial concentration in positive electrode [mol.m-3]"] > 0 ) + + @pytest.mark.unit + def test_evaluate_symbol(self): + parameter_set = pybop.ParameterSet.pybamm("Chen2020") + porosity = parameter_set["Positive electrode porosity"] + assert isinstance(porosity, float) + + for param in [ + 1.0 + porosity, + 1.0 + Scalar(porosity), + 1.0 + Parameter("Positive electrode porosity"), + 1.0 + FunctionParameter("Positive electrode porosity", inputs={}), + ]: + value = pybop.ParameterSet.evaluate_symbol(param, parameter_set) + assert isinstance(value, float) + np.testing.assert_allclose(value, 1.0 + porosity)