From 5955c948ea09165ce8afafbde3a58ac2a5d102d6 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 4 Aug 2023 17:10:49 +0100 Subject: [PATCH 01/18] Progress: Updt. model build structure for parameterisation groups --- examples/Initial_API.py | 22 ++++----- pybop/parameterisation.py | 93 +++++++++++++++++++++++++-------------- 2 files changed, 72 insertions(+), 43 deletions(-) diff --git a/examples/Initial_API.py b/examples/Initial_API.py index f0bb4e4a3..6d37f5157 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -4,21 +4,23 @@ # Form observations Measurements = pd.read_csv("examples/Chen_example.csv", comment='#').to_numpy() -observations = dict( - Time = pybop.Observed(["Time [s]"], Measurements[:,0]), - Current = pybop.Observed(["Current function [A]"], Measurements[:,1]), - Voltage = pybop.Observed(["Voltage [V]"], Measurements[:,2]) -) +observations = { + "Time [s]": pybop.Observed(["Time [s]"], Measurements[:,0]), + "Current function [A]": pybop.Observed(["Current function [A]"], Measurements[:,1]), + "Voltage [V]": pybop.Observed(["Voltage [V]"], Measurements[:,2]) +} # Define model model = pybop.models.lithium_ion.SPM() +model.parameter_values = pybop.ParameterSet("Chen2020") #To implement # Fitting parameters -params = ( - pybop.Parameter("Electrode height [m]", prior = pybop.Gaussian(0,1), bounds = (0.03,0.1)), - pybop.Parameter("Negative particle radius [m]", prior = pybop.Uniform(0,1), bounds = (0.1e-6,0.8e-6)), - pybop.Parameter("Positive particle radius [m]", prior = pybop.Uniform(0,1), bounds = (0.1e-5,0.8e-5)) -) +params = { + "Upper voltage cut-off [V]": pybop.Parameter("Upper voltage cut-off [V]", prior = pybop.Gaussian(4.0,0.01), bounds = [3.8,4.1]) + # "Electrode height [m]": pybop.Parameter("Electrode height [m]", prior = pybop.Gaussian(0,1), bounds = [0.03,0.1]), + # "Negative particle radius [m]": pybop.Parameter("Negative particle radius [m]", prior = pybop.Gaussian(0,1), bounds = [0.1e-6,0.8e-6]), + # "Positive particle radius [m]": pybop.Parameter("Positive particle radius [m]", prior = pybop.Gaussian(0,1), bounds = [0.1e-5,0.8e-5]) +} parameterisation = pybop.Parameterisation(model, observations=observations, fit_parameters=params) diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 46c82bca4..79ae2526a 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -20,32 +20,59 @@ def __init__(self, model, observations, fit_parameters, x0=None, options=None): self.observations = observations self.options = options self.bounds = dict( - lower=[p.bounds[0] for p in fit_parameters], - upper=[p.bounds[1] for p in fit_parameters], + lower=[self.fit_parameters[p].bounds[0] for p in self.fit_parameters], + upper=[self.fit_parameters[p].bounds[1] for p in self.fit_parameters], ) - if options is not None: - self.parameter_set = options.parameter_set + if model.pybamm_model is not None: + if options is not None: + self.parameter_set = options.parameter_set + else: + self.parameter_set = self.model.pybamm_model.default_parameter_values + + try: + self.parameter_set["Current function [A]"] = pybamm.Interpolant( + self.observations["Time [s]"].data, + self.observations["Current function [A]"].data, + pybamm.t, + ) + except: + raise ValueError("Current function not supplied") + + if x0 is None: + self.x0 = np.mean([self.bounds["lower"], self.bounds["upper"]], axis=0) + + # self.sim = pybop.Simulation( + # self.model.pybamm_model, + # parameter_values=self.parameter_set, + # solver=pybamm.CasadiSolver(), + # ) + + # set input parameters in parameter set from fitting parameters + for i in self.fit_parameters: + self.parameter_set[i] = "[input]" + + self.fit_dict = {p: self.fit_parameters[p].prior.mean for p in self.fit_parameters} + + #Set up geometry on model + geometry = self.model.pybamm_model.default_geometry + + # Set up parameters for geometry and model + self.parameter_set.process_model(self.model.pybamm_model) + self.parameter_set.process_geometry(geometry) + + # Mesh the model + mesh = pybamm.Mesh(geometry, self.model.pybamm_model.default_submesh_types, self.model.pybamm_model.default_var_pts) #update + + # Discretise the model + disc = pybamm.Discretisation(mesh, self.model.pybamm_model.default_spatial_methods) + disc.process_model(self.model.pybamm_model) + + # Set solver + self.solver = pybamm.CasadiSolver() else: - self.parameter_set = self.model.pybamm_model.default_parameter_values - - try: - self.parameter_set["Current function [A]"] = pybamm.Interpolant( - self.observations["Time"].data, - self.observations["Current"].data, - pybamm.t, - ) - except: - raise ValueError("Current function not supplied") - - if x0 is None: - self.x0 = np.mean([self.bounds["lower"], self.bounds["upper"]], axis=0) - - self.sim = pybop.Simulation( - self.model.pybamm_model, - parameter_values=self.parameter_set, - solver=pybamm.CasadiSolver(), - ) + raise ValueError("No pybamm model supplied") + def map(self, x0): """ @@ -71,15 +98,16 @@ def rmse(self, method=None): """ def step(x, grad): - for i in range(len(self.fit_parameters)): - self.parameter_set.update({self.fit_parameters[i].name: x[i]}) + for i,p in enumerate(self.fit_dict): + self.fit_dict[p] = x[i] + print(self.fit_dict) - self.sim = pybamm.Simulation( - self.model.pybamm_model, - parameter_values=self.parameter_set, - solver=pybamm.CasadiSolver(), - ) - y_hat = self.sim.solve()["Terminal voltage [V]"].data + # self.sim = pybamm.Simulation( + # self.model.pybamm_model, + # parameter_values=self.parameter_set, + # solver=pybamm.CasadiSolver(), + # ) + y_hat = self.solver.solve(self.model.pybamm_model, inputs=self.fit_dict)["Terminal voltage [V]"].data try: res = np.sqrt( @@ -103,5 +131,4 @@ def mle(self, method=None): """ Maximum likelihood estimation. """ - pass - + pass \ No newline at end of file From fb082fc7bbb9e93f367064c9be26bc5f55887ab6 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Sat, 5 Aug 2023 18:38:04 +0100 Subject: [PATCH 02/18] Progress: Add parameter_sets class with model class inclusion --- examples/Initial_API.py | 2 +- pybop/__init__.py | 4 ++++ pybop/models/lithium_ion/spm.py | 3 +-- pybop/parameter_sets.py | 16 ++++++++++++++++ pybop/parameterisation.py | 4 ++-- 5 files changed, 24 insertions(+), 5 deletions(-) create mode 100644 pybop/parameter_sets.py diff --git a/examples/Initial_API.py b/examples/Initial_API.py index 6d37f5157..48749344f 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -12,7 +12,7 @@ # Define model model = pybop.models.lithium_ion.SPM() -model.parameter_values = pybop.ParameterSet("Chen2020") #To implement +model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") #To implement # Fitting parameters params = { diff --git a/pybop/__init__.py b/pybop/__init__.py index 7cc8690e3..14204e866 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -30,6 +30,10 @@ # from .models import lithium_ion +# +# Parameterisation class +# +from .parameter_sets import ParameterSet # # Parameterisation class # diff --git a/pybop/models/lithium_ion/spm.py b/pybop/models/lithium_ion/spm.py index 367596ddd..99b118a99 100644 --- a/pybop/models/lithium_ion/spm.py +++ b/pybop/models/lithium_ion/spm.py @@ -4,11 +4,10 @@ class SPM: """ - Composition of the SPM class in PyBaMM. - """ def __init__(self): self.pybamm_model = pybamm.lithium_ion.SPM() self.name = "Single Particle Model" + self.parameter_set = self.pybamm_model.default_parameter_values diff --git a/pybop/parameter_sets.py b/pybop/parameter_sets.py new file mode 100644 index 000000000..299b1d41c --- /dev/null +++ b/pybop/parameter_sets.py @@ -0,0 +1,16 @@ +import pybamm +import pybop + +class ParameterSet: + """ + Class for creating parameter sets in pybop. + """ + + def __new__(cls, method, name): + if method.casefold() == "pybamm": + try: + return pybamm.ParameterValues(name).copy() + except: + raise ValueError("Parameter set not found") + else: + raise ValueError("Only PybaMM parameter sets are currently implemented") \ No newline at end of file diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 79ae2526a..ee7672ae6 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -25,8 +25,8 @@ def __init__(self, model, observations, fit_parameters, x0=None, options=None): ) if model.pybamm_model is not None: - if options is not None: - self.parameter_set = options.parameter_set + if model.parameter_set is not None: + self.parameter_set = model.parameter_set else: self.parameter_set = self.model.pybamm_model.default_parameter_values From d332a3f0dec745253dff09e23e3e794eefb3e95b Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 9 Aug 2023 09:02:09 +0100 Subject: [PATCH 03/18] Updt. parameterisation for model_build --- pybop/parameterisation.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index ee7672ae6..5f1f591dd 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -72,7 +72,24 @@ def __init__(self, model, observations, fit_parameters, x0=None, options=None): self.solver = pybamm.CasadiSolver() else: raise ValueError("No pybamm model supplied") - + + def build_model(self): + """ + Build the model. + """ + + if self.model.pybamm_model.built_model: + return + elif self.model.pybamm_model.is_discretised: + self.model.pybamm_model._model_with_set_params = self.model.pybamm_model + self.model.pybamm_model._built_model = self.pybamm_model + else: + # Set parameters + # Mesh + # Discretise + # Built model + pass + def map(self, x0): """ From 52362fd3354f30238bfd64138a50ab62fb396c43 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 9 Aug 2023 18:00:21 +0100 Subject: [PATCH 04/18] Updt parameterisation cls w/ build_model(), set_init_soc(), set_params() --- examples/Initial_API.py | 3 +- pybop/__init__.py | 1 + pybop/parameter_sets.py | 3 +- pybop/parameterisation.py | 173 +++++++++------ pybop/simulation.py | 448 -------------------------------------- 5 files changed, 112 insertions(+), 516 deletions(-) delete mode 100644 pybop/simulation.py diff --git a/examples/Initial_API.py b/examples/Initial_API.py index 48749344f..5e3efee98 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -16,8 +16,9 @@ # Fitting parameters params = { - "Upper voltage cut-off [V]": pybop.Parameter("Upper voltage cut-off [V]", prior = pybop.Gaussian(4.0,0.01), bounds = [3.8,4.1]) + # "Upper voltage cut-off [V]": pybop.Parameter("Upper voltage cut-off [V]", prior = pybop.Gaussian(4.0,0.01), bounds = [3.8,4.1]) # "Electrode height [m]": pybop.Parameter("Electrode height [m]", prior = pybop.Gaussian(0,1), bounds = [0.03,0.1]), + "Negative electrode Bruggeman coefficient (electrolyte)": pybop.Parameter("Negative electrode Bruggeman coefficient (electrolyte)", prior = pybop.Gaussian(1.5,0.1), bounds = [0.8,1.7]), # "Negative particle radius [m]": pybop.Parameter("Negative particle radius [m]", prior = pybop.Gaussian(0,1), bounds = [0.1e-6,0.8e-6]), # "Positive particle radius [m]": pybop.Parameter("Positive particle radius [m]", prior = pybop.Gaussian(0,1), bounds = [0.1e-5,0.8e-5]) } diff --git a/pybop/__init__.py b/pybop/__init__.py index 14204e866..b32cce831 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -34,6 +34,7 @@ # Parameterisation class # from .parameter_sets import ParameterSet + # # Parameterisation class # diff --git a/pybop/parameter_sets.py b/pybop/parameter_sets.py index 299b1d41c..eb84b9da1 100644 --- a/pybop/parameter_sets.py +++ b/pybop/parameter_sets.py @@ -1,6 +1,7 @@ import pybamm import pybop + class ParameterSet: """ Class for creating parameter sets in pybop. @@ -13,4 +14,4 @@ def __new__(cls, method, name): except: raise ValueError("Parameter set not found") else: - raise ValueError("Only PybaMM parameter sets are currently implemented") \ No newline at end of file + raise ValueError("Only PybaMM parameter sets are currently implemented") diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 5f1f591dd..6dcf6b89d 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -14,82 +14,125 @@ class Parameterisation: Parameterisation class for pybop. """ - def __init__(self, model, observations, fit_parameters, x0=None, options=None): - self.model = model + def __init__( + self, + model, + observations, + fit_parameters, + geometry=None, + submesh_types=None, + var_pts=None, + spatial_methods=None, + solver=None, + x0=None, + check_model=True, + init_soc=None, + ): + self.model = model.pybamm_model + self._unprocessed_model = model.pybamm_model self.fit_parameters = fit_parameters self.observations = observations - self.options = options self.bounds = dict( lower=[self.fit_parameters[p].bounds[0] for p in self.fit_parameters], upper=[self.fit_parameters[p].bounds[1] for p in self.fit_parameters], ) + self._model_with_set_params = None + self._built_model = None + self._geometry = geometry or self.model.default_geometry + self._submesh_types = submesh_types or self.model.default_submesh_types + self._var_pts = var_pts or self.model.default_var_pts + self._spatial_methods = spatial_methods or self.model.default_spatial_methods + self.solver = solver or self.model.default_solver + + if model.parameter_set is not None: + self.parameter_set = model.parameter_set + else: + self.parameter_set = self.model.default_parameter_values - if model.pybamm_model is not None: - if model.parameter_set is not None: - self.parameter_set = model.parameter_set - else: - self.parameter_set = self.model.pybamm_model.default_parameter_values - - try: - self.parameter_set["Current function [A]"] = pybamm.Interpolant( - self.observations["Time [s]"].data, - self.observations["Current function [A]"].data, - pybamm.t, - ) - except: - raise ValueError("Current function not supplied") + try: + self.parameter_set["Current function [A]"] = pybamm.Interpolant( + self.observations["Time [s]"].data, + self.observations["Current function [A]"].data, + pybamm.t, + ) + except: + raise ValueError("Current function not supplied") - if x0 is None: - self.x0 = np.mean([self.bounds["lower"], self.bounds["upper"]], axis=0) + if x0 is None: + self.x0 = np.mean([self.bounds["lower"], self.bounds["upper"]], axis=0) - # self.sim = pybop.Simulation( - # self.model.pybamm_model, - # parameter_values=self.parameter_set, - # solver=pybamm.CasadiSolver(), - # ) + # set input parameters in parameter set from fitting parameters + for i in self.fit_parameters: + self.parameter_set[i] = "[input]" - # set input parameters in parameter set from fitting parameters - for i in self.fit_parameters: - self.parameter_set[i] = "[input]" - - self.fit_dict = {p: self.fit_parameters[p].prior.mean for p in self.fit_parameters} + self._unprocessed_parameter_set = self.parameter_set + self.fit_dict = { + p: self.fit_parameters[p].prior.mean for p in self.fit_parameters + } - #Set up geometry on model - geometry = self.model.pybamm_model.default_geometry + # Set t_eval + self.time_data = self.parameter_set["Current function [A]"].x[0] - # Set up parameters for geometry and model - self.parameter_set.process_model(self.model.pybamm_model) - self.parameter_set.process_geometry(geometry) + self.build_model(check_model=check_model, init_soc=init_soc) - # Mesh the model - mesh = pybamm.Mesh(geometry, self.model.pybamm_model.default_submesh_types, self.model.pybamm_model.default_var_pts) #update + # Set solver + self.solver = pybamm.CasadiSolver() - # Discretise the model - disc = pybamm.Discretisation(mesh, self.model.pybamm_model.default_spatial_methods) - disc.process_model(self.model.pybamm_model) + def build_model(self, check_model=True, init_soc=None): + """ + Build the model (if not built already). + """ + if init_soc is not None: + self.set_init_soc(init_soc) # define this function - # Set solver - self.solver = pybamm.CasadiSolver() + if self._built_model: + return + elif self.model.is_discretised: + self.model._model_with_set_params = self.model.pybamm_model + self.model._built_model = self.pybamm_model else: - raise ValueError("No pybamm model supplied") - - def build_model(self): + self.set_params() + self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) + self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods) + self._built_model = self._disc.process_model( + self._model_with_set_params, inplace=False, check_model=check_model + ) + + # Clear solver + self.solver._model_set_up = {} + + def set_init_soc(self, init_soc): + """ + Set the initial state of charge. + """ + if self._built_initial_soc != init_soc: + # reset + self._model_with_set_params = None + self._built_model = None + self.op_conds_to_built_models = None + self.op_conds_to_built_solvers = None + + param = self.model.param + self.parameter_set = ( + self._unprocessed_parameter_set.set_initial_stoichiometries( + init_soc, param=param, inplace=False + ) + ) + # Save solved initial SOC in case we need to re-build the model + self._built_initial_soc = init_soc + + def set_params(self): """ - Build the model. + Set the parameters in the model. """ - - if self.model.pybamm_model.built_model: + if self._model_with_set_params: return - elif self.model.pybamm_model.is_discretised: - self.model.pybamm_model._model_with_set_params = self.model.pybamm_model - self.model.pybamm_model._built_model = self.pybamm_model - else: - # Set parameters - # Mesh - # Discretise - # Built model - pass - + + self._model_with_set_params = self.parameter_set.process_model( + self._unprocessed_model, inplace=False + ) + self.parameter_set.process_geometry(self._geometry) + self.model = self._model_with_set_params def map(self, x0): """ @@ -115,25 +158,23 @@ def rmse(self, method=None): """ def step(x, grad): - for i,p in enumerate(self.fit_dict): + for i, p in enumerate(self.fit_dict): self.fit_dict[p] = x[i] print(self.fit_dict) - # self.sim = pybamm.Simulation( - # self.model.pybamm_model, - # parameter_values=self.parameter_set, - # solver=pybamm.CasadiSolver(), - # ) - y_hat = self.solver.solve(self.model.pybamm_model, inputs=self.fit_dict)["Terminal voltage [V]"].data + y_hat = self.solver.solve( + self._built_model, self.time_data, inputs=self.fit_dict + )["Terminal voltage [V]"].data try: res = np.sqrt( - np.mean((self.observations["Voltage"].data[1] - y_hat) ** 2) + np.mean((self.observations["Voltage [V]"].data[1] - y_hat) ** 2) ) except: raise ValueError( "Measurement and modelled data length mismatch, potentially due to voltage cut-offs" ) + print(res) return res if method == "nlopt": @@ -148,4 +189,4 @@ def mle(self, method=None): """ Maximum likelihood estimation. """ - pass \ No newline at end of file + pass diff --git a/pybop/simulation.py b/pybop/simulation.py deleted file mode 100644 index c1c33b9b0..000000000 --- a/pybop/simulation.py +++ /dev/null @@ -1,448 +0,0 @@ -import pybamm -import numpy as np -import copy -import pickle -import sys -from functools import lru_cache -import warnings - - -def is_notebook(): - try: - shell = get_ipython().__class__.__name__ - if shell == "ZMQInteractiveShell": # pragma: no cover - # Jupyter notebook or qtconsole - cfg = get_ipython().config - nb = len(cfg["InteractiveShell"].keys()) == 0 - return nb - elif shell == "TerminalInteractiveShell": # pragma: no cover - return False # Terminal running IPython - elif shell == "Shell": # pragma: no cover - return True # Google Colab notebook - else: # pragma: no cover - return False # Other type (?) - except NameError: - return False # Probably standard Python interpreter - - -class Simulation: - """ - - This class constructs the PyBOP simulation class. It was initially built from the PyBaMM simulation class. - - Parameters: - ================ - - - """ - - def __init__( - self, - model, - measured_experiment=None, - geometry=None, - parameter_values=None, - submesh_types=None, - var_pts=None, - spatial_methods=None, - solver=None, - output_variables=None, - C_rate=None, - ): - self.parameter_values = parameter_values or model.default_parameter_values - - # Check to see that current is provided as a drive_cycle - current = self._parameter_values.get("Current function [A]") - if isinstance(current, pybamm.Interpolant): - self.operating_mode = "drive cycle" - elif isinstance(current, pybamm.Interpolant): - # This requires syncing the sampling frequency to ensure equivalent vector lengths - self.operating_mode = "without experiment" - if C_rate: - self.C_rate = C_rate - self._parameter_values.update( - { - "Current function [A]": self.C_rate - * self._parameter_values["Nominal cell capacity [A.h]"] - } - ) - else: - raise TypeError( - "measured_experiment must be drive_cycle or C_rate with" - "matching sampling frequency between t_eval and measured data" - ) - - self._unprocessed_model = model - self.model = model - self.geometry = geometry or self.model.default_geometry - self.submesh_types = submesh_types or self.model.default_submesh_types - self.var_pts = var_pts or self.model.default_var_pts - self.spatial_methods = spatial_methods or self.model.default_spatial_methods - self.solver = solver or self.model.default_solver - self.output_variables = output_variables - - # Initialize empty built states - self._model_with_set_params = None - self._built_model = None - self._built_initial_soc = None - self.op_conds_to_built_models = None - self.op_conds_to_built_solvers = None - self._mesh = None - self._disc = None - self._solution = None - self.quick_plot = None - - # ignore runtime warnings in notebooks - if is_notebook(): # pragma: no cover - import warnings - - warnings.filterwarnings("ignore") - - self.get_esoh_solver = lru_cache()(self._get_esoh_solver) - - def set_parameters(self): - """ - Setter for parameter values - - Inputs: - ============ - param: The parameter object to set - """ - if self.model_with_set_params: - return - - self._model_with_set_params = self._parameter_values.process_model( - self._unprocessed_model, inplace=False - ) - self._parameter_values.process_geometry(self.geometry) - self.model = self._model_with_set_params - - def build(self, check_model=True, initial_soc=None): - """ - A method to build the model into a system of matrices and vectors suitable for - performing numerical computations. If the model has already been built or - solved then this function will have no effect. This method will automatically set the parameters - if they have not already been set. - - Parameters - ---------- - check_model : bool, optional - If True, model checks are performed after discretisation (see - :meth:`pybamm.Discretisation.process_model`). Default is True. - initial_soc : float, optional - Initial State of Charge (SOC) for the simulation. Must be between 0 and 1. - If given, overwrites the initial concentrations provided in the parameter - set. - """ - if initial_soc is not None: - self.set_initial_soc(initial_soc) - - if self.built_model: - return - elif self.model.is_discretised: - self._model_with_set_params = self.model - self._built_model = self.model - else: - self.set_parameters() - self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) - self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods) - self._built_model = self._disc.process_model( - self._model_with_set_params, inplace=False, check_model=check_model - ) - # rebuilt model so clear solver setup - self._solver._model_set_up = {} - - def setup_for_parameterisation(): - """ - A method to setup self.model for the parameterisation experiment - """ - - def plot(self, output_variables=None, **kwargs): - """ - A method to quickly plot the outputs of the simulation. Creates a - :class:`pybamm.QuickPlot` object (with keyword arguments 'kwargs') and - then calls :meth:`pybamm.QuickPlot.dynamic_plot`. - - Parameters - ---------- - output_variables: list, optional - A list of the variables to plot. - **kwargs - Additional keyword arguments passed to - :meth:`pybamm.QuickPlot.dynamic_plot`. - For a list of all possible keyword arguments see :class:`pybamm.QuickPlot`. - """ - - if self._solution is None: - raise ValueError( - "Model has not been solved, please solve the model before plotting." - ) - - if output_variables is None: - output_variables = self.output_variables - - self.quick_plot = pybop.dynamic_plot( - self._solution, output_variables=output_variables, **kwargs - ) - - return self.quick_plot - - def create_gif(self, number_of_images=80, duration=0.1, output_filename="plot.gif"): - """ - Create a gif of the parameterisation steps created by :meth:`pybamm.Simulation.plot`. - - Parameters - ---------- - number_of_images : int (optional) - Number of images/plots to be compiled for a GIF. - duration : float (optional) - Duration of visibility of a single image/plot in the created GIF. - output_filename : str (optional) - Name of the generated GIF file. - - """ - if self.quick_plot is None: - self.quick_plot = pybamm.QuickPlot(self._solution) - - # create_git needs to be updated - self.quick_plot.create_gif( - number_of_images=number_of_images, - duration=duration, - output_filename=output_filename, - ) - - def solve( - self, - t_eval=None, - solver=None, - check_model=True, - calc_esoh=True, - starting_solution=None, - initial_soc=None, - callbacks=None, - showprogress=False, - **kwargs, - ): - """ - A method to solve the model for parameterisation. This method will automatically build - and set the model parameters if not already done so. - - Parameters - ---------- - t_eval : numeric type, optional - The times (in seconds) at which to compute the solution. Can be - provided as an array of times at which to return the solution, or as a - list `[t0, tf]` where `t0` is the initial time and `tf` is the final time. - If provided as a list the solution is returned at 100 points within the - interval `[t0, tf]`. - - If None and the parameter "Current function [A]" is read from data - (i.e. drive cycle simulation) the model will be solved at the times - provided in the data. - solver : :class:`pybop.BaseSolver`, optional - The solver to use to solve the model. If None, Simulation.solver is used - check_model : bool, optional - If True, model checks are performed after discretisation (see - :meth:`pybamm.Discretisation.process_model`). Default is True. - calc_esoh : bool, optional - Whether to include eSOH variables in the summary variables. If `False` - then only summary variables that do not require the eSOH calculation - are calculated. Default is True. - starting_solution : :class:`pybamm.Solution` - The solution to start stepping from. If None (default), then self._solution - is used. Must be None if not using an experiment. - initial_soc : float, optional - Initial State of Charge (SOC) for the simulation. Must be between 0 and 1. - If given, overwrites the initial concentrations provided in the parameter - set. - callbacks : list of callbacks, optional - A list of callbacks to be called at each time step. Each callback must - implement all the methods defined in :class:`pybamm.callbacks.BaseCallback`. - showprogress : bool, optional - Whether to show a progress bar for cycling. If true, shows a progress bar - for cycles. Has no effect when not used with an experiment. - Default is False. - **kwargs - Additional key-word arguments passed to `solver.solve`. - See :meth:`pybamm.BaseSolver.solve`. - """ - # Setup - if solver is None: - solver = self.solver - - callbacks = pybamm.callbacks.setup_callbacks(callbacks) - logs = {} - - self.build(check_model=check_model, initial_soc=initial_soc) - - if self.operating_mode == "drive cycle": - # For drive cycles (current provided as data) we perform additional - # tests on t_eval (if provided) to ensure the returned solution - # captures the input. - time_data = self._parameter_values["Current function [A]"].x[0] - # If no t_eval is provided, we use the times provided in the data. - if t_eval is None: - pybamm.logger.info("Setting t_eval as specified by the data") - t_eval = time_data - # If t_eval is provided we first check if it contains all of the - # times in the data to within 10-12. If it doesn't, we then check - # that the largest gap in t_eval is smaller than the smallest gap in - # the time data (to ensure the resolution of t_eval is fine enough). - # We only raise a warning here as users may genuinely only want - # the solution returned at some specified points. - elif ( - set(np.round(time_data, 12)).issubset(set(np.round(t_eval, 12))) - ) is False: - warnings.warn( - """ - t_eval does not contain all of the time points in the data - set. Note: passing t_eval = None automatically sets t_eval - to be the points in the data. - """, - pybamm.SolverWarning, - ) - dt_data_min = np.min(np.diff(time_data)) - dt_eval_max = np.max(np.diff(t_eval)) - if dt_eval_max > dt_data_min + sys.float_info.epsilon: - warnings.warn( - """ - The largest timestep in t_eval ({}) is larger than - the smallest timestep in the data ({}). The returned - solution may not have the correct resolution to accurately - capture the input. Try refining t_eval. Alternatively, - passing t_eval = None automatically sets t_eval to be the - points in the data. - """.format( - dt_eval_max, dt_data_min - ), - pybamm.SolverWarning, - ) - - self._solution = solver.solve(self.built_model, t_eval, **kwargs) - - return self.solution - - def _get_esoh_solver(self, calc_esoh): - if ( - calc_esoh is False - or isinstance(self.model, pybamm.lead_acid.BaseModel) - or isinstance(self.model, pybamm.equivalent_circuit.Thevenin) - or self.model.options["working electrode"] != "both" - ): - return None - - return pybamm.lithium_ion.ElectrodeSOHSolver( - self.parameter_values, self.model.param - ) - - def save(self, filename): - """Save simulation using pickle""" - if self.model.convert_to_format == "python": - # We currently cannot save models in the 'python' format - raise NotImplementedError( - """ - Cannot save simulation if model format is python. - Set model.convert_to_format = 'casadi' instead. - """ - ) - # Clear solver problem (not pickle-able, will automatically be recomputed) - if ( - isinstance(self._solver, pybamm.CasadiSolver) - and self._solver.integrator_specs != {} - ): - self._solver.integrator_specs = {} - - if self.op_conds_to_built_solvers is not None: - for solver in self.op_conds_to_built_solvers.values(): - if ( - isinstance(solver, pybamm.CasadiSolver) - and solver.integrator_specs != {} - ): - solver.integrator_specs = {} - - with open(filename, "wb") as f: - pickle.dump(self, f, pickle.HIGHEST_PROTOCOL) - - def load_sim(filename): - """Load a saved simulation""" - return pybamm.load(filename) - - @property - def model(self): - return self._model - - @model.setter - def model(self, model): - self._model = copy.copy(model) - - @property - def model_with_set_params(self): - return self._model_with_set_params - - @property - def built_model(self): - return self._built_model - - @property - def geometry(self): - return self._geometry - - @geometry.setter - def geometry(self, geometry): - self._geometry = geometry.copy() - - @property - def parameter_values(self): - return self._parameter_values - - @parameter_values.setter - def parameter_values(self, parameter_values): - self._parameter_values = parameter_values.copy() - - @property - def submesh_types(self): - return self._submesh_types - - @submesh_types.setter - def submesh_types(self, submesh_types): - self._submesh_types = submesh_types.copy() - - @property - def mesh(self): - return self._mesh - - @property - def var_pts(self): - return self._var_pts - - @var_pts.setter - def var_pts(self, var_pts): - self._var_pts = var_pts.copy() - - @property - def spatial_methods(self): - return self._spatial_methods - - @spatial_methods.setter - def spatial_methods(self, spatial_methods): - self._spatial_methods = spatial_methods.copy() - - @property - def solver(self): - return self._solver - - @solver.setter - def solver(self, solver): - self._solver = solver.copy() - - @property - def output_variables(self): - return self._output_variables - - @output_variables.setter - def output_variables(self, output_variables): - self._output_variables = copy.copy(output_variables) - - @property - def solution(self): - return self._solution From 57d302558ee2c8de49f9fb8d56066fe1b790c04c Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 10 Aug 2023 17:02:47 +0100 Subject: [PATCH 05/18] Non-geometric parameterisation fitting complete, rm simulation class from init --- examples/Initial_API.py | 9 +++------ pybop/__init__.py | 5 ----- pybop/parameterisation.py | 12 ++++++------ 3 files changed, 9 insertions(+), 17 deletions(-) diff --git a/examples/Initial_API.py b/examples/Initial_API.py index 5e3efee98..1e5ffb7d9 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -12,15 +12,12 @@ # Define model model = pybop.models.lithium_ion.SPM() -model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") #To implement +model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") # Fitting parameters params = { - # "Upper voltage cut-off [V]": pybop.Parameter("Upper voltage cut-off [V]", prior = pybop.Gaussian(4.0,0.01), bounds = [3.8,4.1]) - # "Electrode height [m]": pybop.Parameter("Electrode height [m]", prior = pybop.Gaussian(0,1), bounds = [0.03,0.1]), - "Negative electrode Bruggeman coefficient (electrolyte)": pybop.Parameter("Negative electrode Bruggeman coefficient (electrolyte)", prior = pybop.Gaussian(1.5,0.1), bounds = [0.8,1.7]), - # "Negative particle radius [m]": pybop.Parameter("Negative particle radius [m]", prior = pybop.Gaussian(0,1), bounds = [0.1e-6,0.8e-6]), - # "Positive particle radius [m]": pybop.Parameter("Positive particle radius [m]", prior = pybop.Gaussian(0,1), bounds = [0.1e-5,0.8e-5]) + "Negative electrode active material volume fraction": pybop.Parameter("Negative electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.1,0.9]), + "Positive electrode active material volume fraction": pybop.Parameter("Positive electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.1,0.9]), } parameterisation = pybop.Parameterisation(model, observations=observations, fit_parameters=params) diff --git a/pybop/__init__.py b/pybop/__init__.py index b32cce831..e154c23ff 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -52,11 +52,6 @@ # from .priors import Gaussian, Uniform, Exponential -# -# Simulation class -# -from .simulation import Simulation - # # Optimisation class # diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 6dcf6b89d..90d8df8da 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -5,7 +5,6 @@ # To Do: # checks on observations and parameters -# implement method to update parameterisation without reconstructing the simulation # Geometric parameters might always require reconstruction (i.e. electrode height) @@ -44,6 +43,9 @@ def __init__( self._spatial_methods = spatial_methods or self.model.default_spatial_methods self.solver = solver or self.model.default_solver + # Set solver + self.solver = pybamm.CasadiSolver() + if model.parameter_set is not None: self.parameter_set = model.parameter_set else: @@ -75,9 +77,6 @@ def __init__( self.build_model(check_model=check_model, init_soc=init_soc) - # Set solver - self.solver = pybamm.CasadiSolver() - def build_model(self, check_model=True, init_soc=None): """ Build the model (if not built already). @@ -88,8 +87,8 @@ def build_model(self, check_model=True, init_soc=None): if self._built_model: return elif self.model.is_discretised: - self.model._model_with_set_params = self.model.pybamm_model - self.model._built_model = self.pybamm_model + self.model._model_with_set_params = self.model + self.model._built_model = self.model else: self.set_params() self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) @@ -180,6 +179,7 @@ def step(x, grad): if method == "nlopt": optim = pybop.NLoptOptimize(x0=self.x0) results = optim.optimise(step, self.x0, self.bounds) + else: optim = pybop.NLoptOptimize(method) results = optim.optimise(step, self.x0, self.bounds) From 5df369602dcd83a31a8b2361c7ebfe782508d943 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 16 Aug 2023 12:59:38 +0100 Subject: [PATCH 06/18] Add optim step method w/ lambda, Add initial test framework w/ nox + pytest --- .github/workflows/test_on_push.yaml | 26 +++++++++++++++ examples/Initial_API.py | 8 ++--- noxfile.py | 16 ++++++++++ pybop/parameterisation.py | 49 +++++++++++++++-------------- requirements.txt | 5 --- setup.py | 10 +++--- tests/unit/test_parameterisation.py | 45 ++++++++++++++++++++++++++ 7 files changed, 121 insertions(+), 38 deletions(-) create mode 100644 .github/workflows/test_on_push.yaml create mode 100644 noxfile.py delete mode 100644 requirements.txt create mode 100644 tests/unit/test_parameterisation.py diff --git a/.github/workflows/test_on_push.yaml b/.github/workflows/test_on_push.yaml new file mode 100644 index 000000000..eedb94d68 --- /dev/null +++ b/.github/workflows/test_on_push.yaml @@ -0,0 +1,26 @@ +name: Python package + +on: [push] + +jobs: + build: + + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.8", "3.9", "3.10", "3.11"] + + steps: + - uses: actions/checkout@v3 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install --upgrade pip nox + pip install -e . + - name: Test with pytest + run: | + pytest diff --git a/examples/Initial_API.py b/examples/Initial_API.py index 1e5ffb7d9..c30bbf8ed 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -16,14 +16,14 @@ # Fitting parameters params = { - "Negative electrode active material volume fraction": pybop.Parameter("Negative electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.1,0.9]), - "Positive electrode active material volume fraction": pybop.Parameter("Positive electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.1,0.9]), + "Negative electrode active material volume fraction": pybop.Parameter("Negative electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), + "Positive electrode active material volume fraction": pybop.Parameter("Positive electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), } parameterisation = pybop.Parameterisation(model, observations=observations, fit_parameters=params) -# get RMSE estimate -results, last_optim, num_evals = parameterisation.rmse(method="nlopt") +# get RMSE estimate using NLOpt +results, last_optim, num_evals = parameterisation.rmse(signal="Voltage [V]", method="nlopt") # get MAP estimate, starting at a random initial point in parameter space # parameterisation.map(x0=[p.sample() for p in params]) diff --git a/noxfile.py b/noxfile.py new file mode 100644 index 000000000..d96c45832 --- /dev/null +++ b/noxfile.py @@ -0,0 +1,16 @@ +import nox + +# nox options +nox.options.reuse_existing_virtualenvs = True + +# @nox.session +# def lint(session): +# session.install('flake8') +# session.run('flake8', 'example.py') + + +@nox.session +def tests(session): + session.run_always('pip', 'install', '-e', '.') + session.install('pytest') + session.run('pytest') \ No newline at end of file diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 90d8df8da..e54aa13d2 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -133,9 +133,28 @@ def set_params(self): self.parameter_set.process_geometry(self._geometry) self.model = self._model_with_set_params + def step(self, signal, x, grad): + for i, p in enumerate(self.fit_dict): + self.fit_dict[p] = x[i] + print(self.fit_dict) + + y_hat = self.solver.solve( + self._built_model, self.time_data, inputs=self.fit_dict + )[signal].data + + try: + res = np.sqrt( + np.mean((self.observations["Voltage [V]"].data[1] - y_hat) ** 2) + ) + except: + raise ValueError( + "Measurement and modelled data length mismatch, potentially due to reaching a voltage cut-off" + ) + return res + def map(self, x0): """ - Max a posteriori estimation. + Maximum a posteriori estimation. """ pass @@ -151,38 +170,20 @@ def pem(self, method=None): """ pass - def rmse(self, method=None): + def rmse(self, signal, method=None): """ Calculate the root mean squared error. """ - def step(x, grad): - for i, p in enumerate(self.fit_dict): - self.fit_dict[p] = x[i] - print(self.fit_dict) - - y_hat = self.solver.solve( - self._built_model, self.time_data, inputs=self.fit_dict - )["Terminal voltage [V]"].data - - try: - res = np.sqrt( - np.mean((self.observations["Voltage [V]"].data[1] - y_hat) ** 2) - ) - except: - raise ValueError( - "Measurement and modelled data length mismatch, potentially due to voltage cut-offs" - ) - print(res) - return res - if method == "nlopt": optim = pybop.NLoptOptimize(x0=self.x0) - results = optim.optimise(step, self.x0, self.bounds) + results = optim.optimise( + lambda x, grad: self.step(signal, x, grad), self.x0, self.bounds + ) else: optim = pybop.NLoptOptimize(method) - results = optim.optimise(step, self.x0, self.bounds) + results = optim.optimise(self.step, self.x0, self.bounds) return results def mle(self, method=None): diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index dae82630d..000000000 --- a/requirements.txt +++ /dev/null @@ -1,5 +0,0 @@ -numpy -pandas -scipy -pybamm -matplotlib \ No newline at end of file diff --git a/setup.py b/setup.py index 42d1d101f..bf5523cde 100644 --- a/setup.py +++ b/setup.py @@ -31,11 +31,11 @@ url='https://github.com/pybop-team/PyBOP', # List of packages to install with this one install_requires=[ - "pybamm>=23.1" - "numpy>=1.16", - "scipy>=1.3", - "pandas>=0.24", - "casadi>=3.6.0", + "pybamm>=23.1", + "numpy>=1.25", + "scipy>=1.11", + "pandas>=2.0", + "casadi>=3.6", "nlopt>=2.6", ], # https://pypi.org/classifiers/ diff --git a/tests/unit/test_parameterisation.py b/tests/unit/test_parameterisation.py new file mode 100644 index 000000000..551366d31 --- /dev/null +++ b/tests/unit/test_parameterisation.py @@ -0,0 +1,45 @@ +import pybop +import pytest +import numpy as np +import pandas as pd + + +class TestParameterisation: + def test_rmse(self): + # Form observations + Measurements = pd.read_csv("examples/Chen_example.csv", comment='#').to_numpy() + observations = { + "Time [s]": pybop.Observed(["Time [s]"], Measurements[:,0]), + "Current function [A]": pybop.Observed(["Current function [A]"], Measurements[:,1]), + "Voltage [V]": pybop.Observed(["Voltage [V]"], Measurements[:,2]) + } + + # Define model + model = pybop.models.lithium_ion.SPM() + model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") + + # Fitting parameters + params = { + "Negative electrode active material volume fraction": pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.1), + bounds=[0.05, 0.95], + ), + "Positive electrode active material volume fraction": pybop.Parameter( + "Positive electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.1), + bounds=[0.05, 0.95], + ), + } + + parameterisation = pybop.Parameterisation( + model, observations=observations, fit_parameters=params + ) + + # get RMSE estimate using NLOpt + results, last_optim, num_evals = parameterisation.rmse( + signal="Voltage [V]", method="nlopt" + ) + # Assertions + np.testing.assert_allclose(last_optim, 1e-1, atol=5e-2) + \ No newline at end of file From 66d58565b228615bfced9d4bfe896de3706aebbe Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 16 Aug 2023 13:09:08 +0100 Subject: [PATCH 07/18] Fix test_on_push for nox --- .github/workflows/test_on_push.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test_on_push.yaml b/.github/workflows/test_on_push.yaml index eedb94d68..9fc78541f 100644 --- a/.github/workflows/test_on_push.yaml +++ b/.github/workflows/test_on_push.yaml @@ -1,4 +1,4 @@ -name: Python package +name: PyBOP on: [push] @@ -21,6 +21,6 @@ jobs: python -m pip install --upgrade pip pip install --upgrade pip nox pip install -e . - - name: Test with pytest + - name: Test with nox run: | - pytest + nox From c0001e5c55886e2ea1b0ba38d1a58e3160921b4c Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 16 Aug 2023 13:12:14 +0100 Subject: [PATCH 08/18] Updt numpy version requirement, add fail-fast false condition --- .github/workflows/test_on_push.yaml | 1 + setup.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test_on_push.yaml b/.github/workflows/test_on_push.yaml index 9fc78541f..04f4afc26 100644 --- a/.github/workflows/test_on_push.yaml +++ b/.github/workflows/test_on_push.yaml @@ -7,6 +7,7 @@ jobs: runs-on: ubuntu-latest strategy: + fail-fast: false matrix: python-version: ["3.8", "3.9", "3.10", "3.11"] diff --git a/setup.py b/setup.py index bf5523cde..67ed0aef4 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ # List of packages to install with this one install_requires=[ "pybamm>=23.1", - "numpy>=1.25", + "numpy>=1.16", "scipy>=1.11", "pandas>=2.0", "casadi>=3.6", From 38bddb072ca36469d3ee2d04abfde4b0f90bf008 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 16 Aug 2023 13:24:03 +0100 Subject: [PATCH 09/18] Updt dependancy requirements for python 3.8 support --- setup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 67ed0aef4..a9a328598 100644 --- a/setup.py +++ b/setup.py @@ -33,8 +33,8 @@ install_requires=[ "pybamm>=23.1", "numpy>=1.16", - "scipy>=1.11", - "pandas>=2.0", + "scipy>=1.3", + "pandas>=1.0", "casadi>=3.6", "nlopt>=2.6", ], From 1fb9fac65ed5d255f36ab3b434020315ec079416 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 16 Aug 2023 15:32:41 +0100 Subject: [PATCH 10/18] Updt README w/ build status --- README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.md b/README.md index da868b90f..441e2c469 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,11 @@ # PyBOP - A *Py*thon package for *B*attery *O*ptimisation and *P*arameterisation +
+ +[![Build Status](https://github.com/pybop-team/PyBOP/actions/workflows/CI/badge.svg)](https://github.com/pybop-team/PyBOP/actions) + +
+ PyBOP aims to be a modular library for the parameterisation and optimisation of battery models, with a particular focus on classes built around [PyBaMM](https://github.com/pybamm-team/PyBaMM) models. The figure below gives the current conceptual idea of PyBOP's structure. This will likely evolve as development progresses.

From cc10421d4dd99209592d70c91f15dab83c396a2b Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Wed, 16 Aug 2023 15:41:07 +0100 Subject: [PATCH 11/18] Test icon syntax fix --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 441e2c469..20f03d5b7 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@

-[![Build Status](https://github.com/pybop-team/PyBOP/actions/workflows/CI/badge.svg)](https://github.com/pybop-team/PyBOP/actions) +[![Build Status](https://github.com/pybop-team/PyBOP/actions/workflows/test_on_push.yaml/badge.svg?branch=develop)](https://github.com/pybop-team/PyBOP/actions/workflows/test_on_push.yaml)
From 6b7ab48a391f3116cedea3c0c7107cb3ac5071a0 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 17 Aug 2023 16:19:40 +0100 Subject: [PATCH 12/18] Updt. x0 to be sampled from priors, Progress: API to list formation --- examples/Initial_API.py | 10 +++++----- pybop/parameterisation.py | 12 ++++++++++-- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/examples/Initial_API.py b/examples/Initial_API.py index c30bbf8ed..261e296ef 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -4,11 +4,11 @@ # Form observations Measurements = pd.read_csv("examples/Chen_example.csv", comment='#').to_numpy() -observations = { - "Time [s]": pybop.Observed(["Time [s]"], Measurements[:,0]), - "Current function [A]": pybop.Observed(["Current function [A]"], Measurements[:,1]), - "Voltage [V]": pybop.Observed(["Voltage [V]"], Measurements[:,2]) -} +observations = [ + pybop.Observed(["Time [s]"], Measurements[:,0]), + pybop.Observed(["Current function [A]"], Measurements[:,1]), + pybop.Observed(["Voltage [V]"], Measurements[:,2]) +] # Define model model = pybop.models.lithium_ion.SPM() diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index e54aa13d2..a31cc31ea 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -30,7 +30,12 @@ def __init__( self.model = model.pybamm_model self._unprocessed_model = model.pybamm_model self.fit_parameters = fit_parameters - self.observations = observations + + self.observations = {o.name: o for o in observations} # needs to be fixed + for name in ["Time [s]", "Current function [A]"]: + if name not in self.observations: + raise ValueError(f"expected {name} in list of observations") + self.bounds = dict( lower=[self.fit_parameters[p].bounds[0] for p in self.fit_parameters], upper=[self.fit_parameters[p].bounds[1] for p in self.fit_parameters], @@ -61,7 +66,10 @@ def __init__( raise ValueError("Current function not supplied") if x0 is None: - self.x0 = np.mean([self.bounds["lower"], self.bounds["upper"]], axis=0) + self.x0 = np.zeros(len(self.fit_parameters)) + + for i, j in enumerate(self.fit_parameters): + self.x0[i] = self.fit_parameters[j].prior.rvs(1) # set input parameters in parameter set from fitting parameters for i in self.fit_parameters: From d50702e284e92cc543c693390e6c98229f76b746 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 17 Aug 2023 18:07:16 +0100 Subject: [PATCH 13/18] API change to observation & params definition to type::list --- examples/Initial_API.py | 14 +++++++------- pybop/parameterisation.py | 6 +++--- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/Initial_API.py b/examples/Initial_API.py index 261e296ef..a884409bb 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -5,9 +5,9 @@ # Form observations Measurements = pd.read_csv("examples/Chen_example.csv", comment='#').to_numpy() observations = [ - pybop.Observed(["Time [s]"], Measurements[:,0]), - pybop.Observed(["Current function [A]"], Measurements[:,1]), - pybop.Observed(["Voltage [V]"], Measurements[:,2]) + pybop.Observed("Time [s]", Measurements[:,0]), + pybop.Observed("Current function [A]", Measurements[:,1]), + pybop.Observed("Voltage [V]", Measurements[:,2]) ] # Define model @@ -15,10 +15,10 @@ model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") # Fitting parameters -params = { - "Negative electrode active material volume fraction": pybop.Parameter("Negative electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), - "Positive electrode active material volume fraction": pybop.Parameter("Positive electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), -} +params = [ + pybop.Parameter("Negative electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), + pybop.Parameter("Positive electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), +] parameterisation = pybop.Parameterisation(model, observations=observations, fit_parameters=params) diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index a31cc31ea..008e38bfd 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -29,9 +29,10 @@ def __init__( ): self.model = model.pybamm_model self._unprocessed_model = model.pybamm_model - self.fit_parameters = fit_parameters + self.fit_parameters = {o.name: o for o in fit_parameters} + self.observations = {o.name: o for o in observations} - self.observations = {o.name: o for o in observations} # needs to be fixed + # Check that the observations contain time and current for name in ["Time [s]", "Current function [A]"]: if name not in self.observations: raise ValueError(f"expected {name} in list of observations") @@ -67,7 +68,6 @@ def __init__( if x0 is None: self.x0 = np.zeros(len(self.fit_parameters)) - for i, j in enumerate(self.fit_parameters): self.x0[i] = self.fit_parameters[j].prior.rvs(1) From ee7fe83f7727672c9ab12a1d2846123d39c25df4 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 18 Aug 2023 09:03:32 +0100 Subject: [PATCH 14/18] refactor API change for observ. & param class type::list, tests update for change --- examples/Initial_API.py | 34 +++++++++++++++++++---------- pybop/parameterisation.py | 2 +- tests/unit/test_parameterisation.py | 21 +++++++++--------- 3 files changed, 34 insertions(+), 23 deletions(-) diff --git a/examples/Initial_API.py b/examples/Initial_API.py index a884409bb..25faa65b7 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -3,30 +3,42 @@ import numpy as np # Form observations -Measurements = pd.read_csv("examples/Chen_example.csv", comment='#').to_numpy() +Measurements = pd.read_csv("examples/Chen_example.csv", comment="#").to_numpy() observations = [ - pybop.Observed("Time [s]", Measurements[:,0]), - pybop.Observed("Current function [A]", Measurements[:,1]), - pybop.Observed("Voltage [V]", Measurements[:,2]) + pybop.Observed("Time [s]", Measurements[:, 0]), + pybop.Observed("Current function [A]", Measurements[:, 1]), + pybop.Observed("Voltage [V]", Measurements[:, 2]), ] - + # Define model model = pybop.models.lithium_ion.SPM() model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") # Fitting parameters params = [ - pybop.Parameter("Negative electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), - pybop.Parameter("Positive electrode active material volume fraction", prior = pybop.Gaussian(0.6,0.1), bounds = [0.05,0.95]), + pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.1), + bounds=[0.05, 0.95], + ), + pybop.Parameter( + "Positive electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.1), + bounds=[0.05, 0.95], + ), ] -parameterisation = pybop.Parameterisation(model, observations=observations, fit_parameters=params) +parameterisation = pybop.Parameterisation( + model, observations=observations, fit_parameters=params +) # get RMSE estimate using NLOpt -results, last_optim, num_evals = parameterisation.rmse(signal="Voltage [V]", method="nlopt") +results, last_optim, num_evals = parameterisation.rmse( + signal="Voltage [V]", method="nlopt" +) # get MAP estimate, starting at a random initial point in parameter space -# parameterisation.map(x0=[p.sample() for p in params]) +# parameterisation.map(x0=[p.sample() for p in params]) # or sample from posterior # parameterisation.sample(1000, n_chains=4, ....) @@ -35,4 +47,4 @@ # parameterisation.sober() -#Optimisation = pybop.optimisation(model, cost=cost, parameters=parameters, observation=observation) \ No newline at end of file +# Optimisation = pybop.optimisation(model, cost=cost, parameters=parameters, observation=observation) diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 008e38bfd..898f87b52 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -69,7 +69,7 @@ def __init__( if x0 is None: self.x0 = np.zeros(len(self.fit_parameters)) for i, j in enumerate(self.fit_parameters): - self.x0[i] = self.fit_parameters[j].prior.rvs(1) + self.x0[i] = self.fit_parameters[j].prior.rvs(1)[0] # set input parameters in parameter set from fitting parameters for i in self.fit_parameters: diff --git a/tests/unit/test_parameterisation.py b/tests/unit/test_parameterisation.py index 551366d31..123363517 100644 --- a/tests/unit/test_parameterisation.py +++ b/tests/unit/test_parameterisation.py @@ -7,30 +7,30 @@ class TestParameterisation: def test_rmse(self): # Form observations - Measurements = pd.read_csv("examples/Chen_example.csv", comment='#').to_numpy() - observations = { - "Time [s]": pybop.Observed(["Time [s]"], Measurements[:,0]), - "Current function [A]": pybop.Observed(["Current function [A]"], Measurements[:,1]), - "Voltage [V]": pybop.Observed(["Voltage [V]"], Measurements[:,2]) - } + Measurements = pd.read_csv("examples/Chen_example.csv", comment="#").to_numpy() + observations = [ + pybop.Observed("Time [s]", Measurements[:, 0]), + pybop.Observed("Current function [A]", Measurements[:, 1]), + pybop.Observed("Voltage [V]", Measurements[:, 2]), + ] # Define model model = pybop.models.lithium_ion.SPM() model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") # Fitting parameters - params = { - "Negative electrode active material volume fraction": pybop.Parameter( + params = [ + pybop.Parameter( "Negative electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.1), bounds=[0.05, 0.95], ), - "Positive electrode active material volume fraction": pybop.Parameter( + pybop.Parameter( "Positive electrode active material volume fraction", prior=pybop.Gaussian(0.6, 0.1), bounds=[0.05, 0.95], ), - } + ] parameterisation = pybop.Parameterisation( model, observations=observations, fit_parameters=params @@ -42,4 +42,3 @@ def test_rmse(self): ) # Assertions np.testing.assert_allclose(last_optim, 1e-1, atol=5e-2) - \ No newline at end of file From d8b8077ffc5f3f735df1a5078c52b4ecbcde039b Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 18 Aug 2023 10:29:17 +0100 Subject: [PATCH 15/18] Refactor parameterisation & spm class, Add BaseModel class --- pybop/__init__.py | 1 + pybop/models/BaseModel.py | 24 +++++++++ pybop/models/lithium_ion/spm.py | 92 ++++++++++++++++++++++++++++++-- pybop/parameterisation.py | 93 +++------------------------------ 4 files changed, 122 insertions(+), 88 deletions(-) create mode 100644 pybop/models/BaseModel.py diff --git a/pybop/__init__.py b/pybop/__init__.py index e154c23ff..7aeec0553 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -28,6 +28,7 @@ # # Model Classes # +from .models import BaseModel from .models import lithium_ion # diff --git a/pybop/models/BaseModel.py b/pybop/models/BaseModel.py new file mode 100644 index 000000000..a526909bf --- /dev/null +++ b/pybop/models/BaseModel.py @@ -0,0 +1,24 @@ +import pybop + + +class BaseModel: + """ + Base class for PyBOP models + """ + + def __init__(self, name="Base Model"): + self.pybamm_model = None + self.name = name + self.parameter_set = None + + def build(self): + """ + Build the model + """ + pass + + def sim(self): + """ + Simulate the model + """ + pass diff --git a/pybop/models/lithium_ion/spm.py b/pybop/models/lithium_ion/spm.py index 99b118a99..5ea73fbc4 100644 --- a/pybop/models/lithium_ion/spm.py +++ b/pybop/models/lithium_ion/spm.py @@ -1,13 +1,99 @@ import pybop import pybamm +from ..BaseModel import BaseModel -class SPM: +class SPM(BaseModel): """ Composition of the SPM class in PyBaMM. """ - def __init__(self): + def __init__( + self, + name="Single Particle Model", + geometry=None, + submesh_types=None, + var_pts=None, + spatial_methods=None, + solver=None, + ): + super().__init__() self.pybamm_model = pybamm.lithium_ion.SPM() - self.name = "Single Particle Model" + self._unprocessed_model = self.pybamm_model + self.name = name self.parameter_set = self.pybamm_model.default_parameter_values + self._model_with_set_params = None + self._built_model = None + self._geometry = geometry or self.pybamm_model.default_geometry + self._submesh_types = submesh_types or self.pybamm_model.default_submesh_types + self._var_pts = var_pts or self.pybamm_model.default_var_pts + self._spatial_methods = ( + spatial_methods or self.pybamm_model.default_spatial_methods + ) + self.solver = solver or self.pybamm_model.default_solver + + def build_model( + self, + check_model=True, + init_soc=None, + ): + """ + Build the model (if not built already). + """ + if init_soc is not None: + self.set_init_soc(init_soc) # define this function + + if self._built_model: + return + + elif self.pybamm_model.is_discretised: + self.pybamm_model._model_with_set_params = self.pybamm_model + self.pybamm_model._built_model = self.pybamm_model + else: + self.set_params() + self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) + self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods) + self._built_model = self._disc.process_model( + self._model_with_set_params, inplace=False, check_model=check_model + ) + + # Clear solver + self.solver._model_set_up = {} + + def set_init_soc(self, init_soc): + """ + Set the initial state of charge. + """ + if self._built_initial_soc != init_soc: + # reset + self._model_with_set_params = None + self._built_model = None + self.op_conds_to_built_models = None + self.op_conds_to_built_solvers = None + + param = self.model.pybamm_model.param + self.parameter_set = ( + self._unprocessed_parameter_set.set_initial_stoichiometries( + init_soc, param=param, inplace=False + ) + ) + # Save solved initial SOC in case we need to re-build the model + self._built_initial_soc = init_soc + + def set_params(self): + """ + Set the parameters in the model. + """ + if self._model_with_set_params: + return + + self._model_with_set_params = self.parameter_set.process_model( + self._unprocessed_model, inplace=False + ) + self.parameter_set.process_geometry(self._geometry) + self.pybamm_model = self._model_with_set_params + + def sim(self): + """ + Simulate the model + """ diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 898f87b52..70cbd7c93 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -3,11 +3,6 @@ import numpy as np -# To Do: -# checks on observations and parameters -# Geometric parameters might always require reconstruction (i.e. electrode height) - - class Parameterisation: """ Parameterisation class for pybop. @@ -18,17 +13,11 @@ def __init__( model, observations, fit_parameters, - geometry=None, - submesh_types=None, - var_pts=None, - spatial_methods=None, - solver=None, x0=None, check_model=True, init_soc=None, ): - self.model = model.pybamm_model - self._unprocessed_model = model.pybamm_model + self.model = model self.fit_parameters = {o.name: o for o in fit_parameters} self.observations = {o.name: o for o in observations} @@ -41,21 +30,11 @@ def __init__( lower=[self.fit_parameters[p].bounds[0] for p in self.fit_parameters], upper=[self.fit_parameters[p].bounds[1] for p in self.fit_parameters], ) - self._model_with_set_params = None - self._built_model = None - self._geometry = geometry or self.model.default_geometry - self._submesh_types = submesh_types or self.model.default_submesh_types - self._var_pts = var_pts or self.model.default_var_pts - self._spatial_methods = spatial_methods or self.model.default_spatial_methods - self.solver = solver or self.model.default_solver - - # Set solver - self.solver = pybamm.CasadiSolver() - - if model.parameter_set is not None: - self.parameter_set = model.parameter_set + + if self.model.parameter_set is not None: + self.parameter_set = self.model.parameter_set else: - self.parameter_set = self.model.default_parameter_values + self.parameter_set = self.model.pybamm_model.default_parameter_values try: self.parameter_set["Current function [A]"] = pybamm.Interpolant( @@ -83,71 +62,15 @@ def __init__( # Set t_eval self.time_data = self.parameter_set["Current function [A]"].x[0] - self.build_model(check_model=check_model, init_soc=init_soc) - - def build_model(self, check_model=True, init_soc=None): - """ - Build the model (if not built already). - """ - if init_soc is not None: - self.set_init_soc(init_soc) # define this function - - if self._built_model: - return - elif self.model.is_discretised: - self.model._model_with_set_params = self.model - self.model._built_model = self.model - else: - self.set_params() - self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) - self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods) - self._built_model = self._disc.process_model( - self._model_with_set_params, inplace=False, check_model=check_model - ) - - # Clear solver - self.solver._model_set_up = {} - - def set_init_soc(self, init_soc): - """ - Set the initial state of charge. - """ - if self._built_initial_soc != init_soc: - # reset - self._model_with_set_params = None - self._built_model = None - self.op_conds_to_built_models = None - self.op_conds_to_built_solvers = None - - param = self.model.param - self.parameter_set = ( - self._unprocessed_parameter_set.set_initial_stoichiometries( - init_soc, param=param, inplace=False - ) - ) - # Save solved initial SOC in case we need to re-build the model - self._built_initial_soc = init_soc - - def set_params(self): - """ - Set the parameters in the model. - """ - if self._model_with_set_params: - return - - self._model_with_set_params = self.parameter_set.process_model( - self._unprocessed_model, inplace=False - ) - self.parameter_set.process_geometry(self._geometry) - self.model = self._model_with_set_params + self.model.build_model(check_model=check_model, init_soc=init_soc) def step(self, signal, x, grad): for i, p in enumerate(self.fit_dict): self.fit_dict[p] = x[i] print(self.fit_dict) - y_hat = self.solver.solve( - self._built_model, self.time_data, inputs=self.fit_dict + y_hat = self.model.solver.solve( + self.model._built_model, self.time_data, inputs=self.fit_dict )[signal].data try: From e3f49c587f604a2633b916f618344b47558bee82 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 18 Aug 2023 11:26:45 +0100 Subject: [PATCH 16/18] Refactor #2 of parameterisation & spm class, Aligned init. parameter conditions for rmse/nlopt --- pybop/models/lithium_ion/spm.py | 28 ++++++++++++++++++-- pybop/parameterisation.py | 47 +++++++++++++-------------------- 2 files changed, 44 insertions(+), 31 deletions(-) diff --git a/pybop/models/lithium_ion/spm.py b/pybop/models/lithium_ion/spm.py index 5ea73fbc4..c42ec0f0e 100644 --- a/pybop/models/lithium_ion/spm.py +++ b/pybop/models/lithium_ion/spm.py @@ -16,12 +16,13 @@ def __init__( var_pts=None, spatial_methods=None, solver=None, + parameter_set=None, ): super().__init__() self.pybamm_model = pybamm.lithium_ion.SPM() self._unprocessed_model = self.pybamm_model self.name = name - self.parameter_set = self.pybamm_model.default_parameter_values + self.parameter_set = parameter_set or self.pybamm_model.default_parameter_values self._model_with_set_params = None self._built_model = None self._geometry = geometry or self.pybamm_model.default_geometry @@ -77,7 +78,7 @@ def set_init_soc(self, init_soc): init_soc, param=param, inplace=False ) ) - # Save solved initial SOC in case we need to re-build the model + # Save solved initial SOC in case we need to rebuild the model self._built_initial_soc = init_soc def set_params(self): @@ -93,6 +94,29 @@ def set_params(self): self.parameter_set.process_geometry(self._geometry) self.pybamm_model = self._model_with_set_params + def build_parameter_set(self, observations, fit_parameters): + """ + Build the parameter set. + """ + + try: + self.parameter_set["Current function [A]"] = pybamm.Interpolant( + observations["Time [s]"].data, + observations["Current function [A]"].data, + pybamm.t, + ) + except: + raise ValueError("Current function not supplied") + + # set input parameters in parameter set from fitting parameters + for i in fit_parameters: + self.parameter_set[i] = "[input]" + + self._unprocessed_parameter_set = self.parameter_set + + # Set t_eval + self.time_data = self.parameter_set["Current function [A]"].x[0] + def sim(self): """ Simulate the model diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 70cbd7c93..3bc9cc27b 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -16,8 +16,10 @@ def __init__( x0=None, check_model=True, init_soc=None, + verbose=False, ): self.model = model + self.fit_dict = {} self.fit_parameters = {o.name: o for o in fit_parameters} self.observations = {o.name: o for o in observations} @@ -26,51 +28,34 @@ def __init__( if name not in self.observations: raise ValueError(f"expected {name} in list of observations") + # Set bounds self.bounds = dict( lower=[self.fit_parameters[p].bounds[0] for p in self.fit_parameters], upper=[self.fit_parameters[p].bounds[1] for p in self.fit_parameters], ) - if self.model.parameter_set is not None: - self.parameter_set = self.model.parameter_set - else: - self.parameter_set = self.model.pybamm_model.default_parameter_values - - try: - self.parameter_set["Current function [A]"] = pybamm.Interpolant( - self.observations["Time [s]"].data, - self.observations["Current function [A]"].data, - pybamm.t, - ) - except: - raise ValueError("Current function not supplied") - + # Sample from prior for x0 if x0 is None: self.x0 = np.zeros(len(self.fit_parameters)) for i, j in enumerate(self.fit_parameters): - self.x0[i] = self.fit_parameters[j].prior.rvs(1)[0] - - # set input parameters in parameter set from fitting parameters - for i in self.fit_parameters: - self.parameter_set[i] = "[input]" - - self._unprocessed_parameter_set = self.parameter_set - self.fit_dict = { - p: self.fit_parameters[p].prior.mean for p in self.fit_parameters - } + self.x0[i] = self.fit_parameters[j].prior.rvs(1)[ + 0 + ] # Updt to capture dimensions per parameter - # Set t_eval - self.time_data = self.parameter_set["Current function [A]"].x[0] + # Align initial guess with sample from prior + for i, j in enumerate(self.fit_parameters): + self.fit_dict[j] = {j: self.x0[i]} + # Build parameter set and model + self.model.build_parameter_set(self.observations, self.fit_parameters) self.model.build_model(check_model=check_model, init_soc=init_soc) - def step(self, signal, x, grad): + def step(self, signal, x, grad, verbose): for i, p in enumerate(self.fit_dict): self.fit_dict[p] = x[i] - print(self.fit_dict) y_hat = self.model.solver.solve( - self.model._built_model, self.time_data, inputs=self.fit_dict + self.model._built_model, self.model.time_data, inputs=self.fit_dict )[signal].data try: @@ -81,6 +66,10 @@ def step(self, signal, x, grad): raise ValueError( "Measurement and modelled data length mismatch, potentially due to reaching a voltage cut-off" ) + + if verbose: + print(self.fit_dict) + return res def map(self, x0): From d21249ba9823ff4535237680ec987a0a60873d96 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 18 Aug 2023 15:38:10 +0100 Subject: [PATCH 17/18] Data generation in test methods, debug parameter initilisation problem for spm+basemodel, parameterisation.py tidy --- examples/Initial_API.py | 12 +-- pybop/models/BaseModel.py | 4 +- pybop/models/lithium_ion/spm.py | 129 +++++++++++++++++++++------- pybop/parameterisation.py | 38 +++++--- tests/unit/test_parameterisation.py | 34 +++++++- 5 files changed, 161 insertions(+), 56 deletions(-) diff --git a/examples/Initial_API.py b/examples/Initial_API.py index 25faa65b7..5a47e27aa 100644 --- a/examples/Initial_API.py +++ b/examples/Initial_API.py @@ -11,20 +11,20 @@ ] # Define model -model = pybop.models.lithium_ion.SPM() -model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") +parameter_set = pybop.ParameterSet("pybamm", "Chen2020") +model = pybop.models.lithium_ion.SPM(parameter_set=parameter_set) # Fitting parameters params = [ pybop.Parameter( "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.1), - bounds=[0.05, 0.95], + prior=pybop.Gaussian(0.75, 0.05), + bounds=[0.65, 0.85], ), pybop.Parameter( "Positive electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.1), - bounds=[0.05, 0.95], + prior=pybop.Gaussian(0.65, 0.05), + bounds=[0.55, 0.75], ), ] diff --git a/pybop/models/BaseModel.py b/pybop/models/BaseModel.py index a526909bf..809a99e28 100644 --- a/pybop/models/BaseModel.py +++ b/pybop/models/BaseModel.py @@ -7,9 +7,9 @@ class BaseModel: """ def __init__(self, name="Base Model"): - self.pybamm_model = None + # self.pybamm_model = None self.name = name - self.parameter_set = None + # self.parameter_set = None def build(self): """ diff --git a/pybop/models/lithium_ion/spm.py b/pybop/models/lithium_ion/spm.py index c42ec0f0e..133a241b6 100644 --- a/pybop/models/lithium_ion/spm.py +++ b/pybop/models/lithium_ion/spm.py @@ -11,30 +11,39 @@ class SPM(BaseModel): def __init__( self, name="Single Particle Model", + parameter_set=None, geometry=None, submesh_types=None, var_pts=None, spatial_methods=None, solver=None, - parameter_set=None, ): super().__init__() self.pybamm_model = pybamm.lithium_ion.SPM() self._unprocessed_model = self.pybamm_model self.name = name + self.parameter_set = parameter_set or self.pybamm_model.default_parameter_values - self._model_with_set_params = None - self._built_model = None - self._geometry = geometry or self.pybamm_model.default_geometry - self._submesh_types = submesh_types or self.pybamm_model.default_submesh_types - self._var_pts = var_pts or self.pybamm_model.default_var_pts - self._spatial_methods = ( + self._unprocessed_parameter_set = self.parameter_set + + self.geometry = geometry or self.pybamm_model.default_geometry + self.submesh_types = submesh_types or self.pybamm_model.default_submesh_types + self.var_pts = var_pts or self.pybamm_model.default_var_pts + self.spatial_methods = ( spatial_methods or self.pybamm_model.default_spatial_methods ) self.solver = solver or self.pybamm_model.default_solver + self._model_with_set_params = None + self._built_model = None + self._built_initial_soc = None + self._mesh = None + self._disc = None + def build_model( self, + observations, + fit_parameters, check_model=True, init_soc=None, ): @@ -42,24 +51,26 @@ def build_model( Build the model (if not built already). """ if init_soc is not None: - self.set_init_soc(init_soc) # define this function + self.set_init_soc(init_soc) if self._built_model: return elif self.pybamm_model.is_discretised: - self.pybamm_model._model_with_set_params = self.pybamm_model - self.pybamm_model._built_model = self.pybamm_model + self._model_with_set_params = self.pybamm_model + self._built_model = self.pybamm_model else: - self.set_params() - self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts) - self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods) + self.set_params(observations, fit_parameters) + self._mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.var_pts) + self._disc = pybamm.Discretisation(self.mesh, self.spatial_methods) self._built_model = self._disc.process_model( self._model_with_set_params, inplace=False, check_model=check_model ) + # Set t_eval + self.time_data = self._parameter_set["Current function [A]"].x[0] # Clear solver - self.solver._model_set_up = {} + self._solver._model_set_up = {} def set_init_soc(self, init_soc): """ @@ -72,7 +83,7 @@ def set_init_soc(self, init_soc): self.op_conds_to_built_models = None self.op_conds_to_built_solvers = None - param = self.model.pybamm_model.param + param = self.pybamm_model.param self.parameter_set = ( self._unprocessed_parameter_set.set_initial_stoichiometries( init_soc, param=param, inplace=False @@ -81,24 +92,13 @@ def set_init_soc(self, init_soc): # Save solved initial SOC in case we need to rebuild the model self._built_initial_soc = init_soc - def set_params(self): + def set_params(self, observations, fit_parameters): """ Set the parameters in the model. """ - if self._model_with_set_params: + if self.model_with_set_params: return - self._model_with_set_params = self.parameter_set.process_model( - self._unprocessed_model, inplace=False - ) - self.parameter_set.process_geometry(self._geometry) - self.pybamm_model = self._model_with_set_params - - def build_parameter_set(self, observations, fit_parameters): - """ - Build the parameter set. - """ - try: self.parameter_set["Current function [A]"] = pybamm.Interpolant( observations["Time [s]"].data, @@ -112,12 +112,77 @@ def build_parameter_set(self, observations, fit_parameters): for i in fit_parameters: self.parameter_set[i] = "[input]" - self._unprocessed_parameter_set = self.parameter_set - - # Set t_eval - self.time_data = self.parameter_set["Current function [A]"].x[0] + self._model_with_set_params = self._parameter_set.process_model( + self._unprocessed_model, inplace=False + ) + self._parameter_set.process_geometry(self.geometry) + self.pybamm_model = self._model_with_set_params def sim(self): """ Simulate the model """ + + @property + def built_model(self): + return self._built_model + + @property + def parameter_set(self): + return self._parameter_set + + @parameter_set.setter + def parameter_set(self, parameter_set): + self._parameter_set = parameter_set.copy() + + @property + def model_with_set_params(self): + return self._model_with_set_params + + @property + def built_model(self): + return self._built_model + + @property + def geometry(self): + return self._geometry + + @geometry.setter + def geometry(self, geometry): + self._geometry = geometry.copy() + + @property + def submesh_types(self): + return self._submesh_types + + @submesh_types.setter + def submesh_types(self, submesh_types): + self._submesh_types = submesh_types.copy() + + @property + def mesh(self): + return self._mesh + + @property + def var_pts(self): + return self._var_pts + + @var_pts.setter + def var_pts(self, var_pts): + self._var_pts = var_pts.copy() + + @property + def spatial_methods(self): + return self._spatial_methods + + @spatial_methods.setter + def spatial_methods(self, spatial_methods): + self._spatial_methods = spatial_methods.copy() + + @property + def solver(self): + return self._solver + + @solver.setter + def solver(self, solver): + self._solver = solver.copy() diff --git a/pybop/parameterisation.py b/pybop/parameterisation.py index 3bc9cc27b..aa6828d3c 100644 --- a/pybop/parameterisation.py +++ b/pybop/parameterisation.py @@ -1,5 +1,4 @@ import pybop -import pybamm import numpy as np @@ -19,6 +18,7 @@ def __init__( verbose=False, ): self.model = model + self.verbose = verbose self.fit_dict = {} self.fit_parameters = {o.name: o for o in fit_parameters} self.observations = {o.name: o for o in observations} @@ -46,28 +46,42 @@ def __init__( for i, j in enumerate(self.fit_parameters): self.fit_dict[j] = {j: self.x0[i]} - # Build parameter set and model - self.model.build_parameter_set(self.observations, self.fit_parameters) - self.model.build_model(check_model=check_model, init_soc=init_soc) + # Build model with observations and fitting_parameters + self.model.build_model( + self.observations, + self.fit_parameters, + check_model=check_model, + init_soc=init_soc, + ) - def step(self, signal, x, grad, verbose): + def step(self, signal, x, grad): for i, p in enumerate(self.fit_dict): self.fit_dict[p] = x[i] y_hat = self.model.solver.solve( - self.model._built_model, self.model.time_data, inputs=self.fit_dict + self.model.built_model, self.model.time_data, inputs=self.fit_dict )[signal].data - try: - res = np.sqrt( - np.mean((self.observations["Voltage [V]"].data[1] - y_hat) ** 2) + print( + "Last Voltage Values:", y_hat[-1], self.observations["Voltage [V]"].data[-1] + ) + + if len(y_hat) != len(self.observations["Voltage [V]"].data): + print( + "len of vectors:", + len(y_hat), + len(self.observations["Voltage [V]"].data), ) - except: raise ValueError( - "Measurement and modelled data length mismatch, potentially due to reaching a voltage cut-off" + "Measurement and simulated data length mismatch, potentially due to reaching a voltage cut-off" ) - if verbose: + try: + res = np.sqrt(np.mean((self.observations["Voltage [V]"].data - y_hat) ** 2)) + except: + print("Error in RMSE calculation") + + if self.verbose: print(self.fit_dict) return res diff --git a/tests/unit/test_parameterisation.py b/tests/unit/test_parameterisation.py index 123363517..a55dae54b 100644 --- a/tests/unit/test_parameterisation.py +++ b/tests/unit/test_parameterisation.py @@ -1,17 +1,43 @@ import pybop import pytest +import pybamm import numpy as np import pandas as pd class TestParameterisation: + def getdata(self, x1, x2): + model = pybamm.lithium_ion.SPM() + params = pybamm.ParameterValues("Chen2020") + + params.update( + { + "Negative electrode active material volume fraction": x1, + "Positive electrode active material volume fraction": x2, + } + ) + experiment = pybamm.Experiment( + [ + ( + "Discharge at 0.5C for 5 minutes (1 second period)", + "Rest for 2 minutes (1 second period)", + "Charge at 0.5C for 2.5 minutes (1 second period)", + "Rest for 2 minutes (1 second period)", + ), + ] + * 2 + ) + sim = pybamm.Simulation(model, experiment=experiment, parameter_values=params) + return sim.solve(initial_soc=0.5) + def test_rmse(self): # Form observations - Measurements = pd.read_csv("examples/Chen_example.csv", comment="#").to_numpy() + solution = self.getdata(0.6, 0.6) + observations = [ - pybop.Observed("Time [s]", Measurements[:, 0]), - pybop.Observed("Current function [A]", Measurements[:, 1]), - pybop.Observed("Voltage [V]", Measurements[:, 2]), + pybop.Observed("Time [s]", solution["Time [s]"].data), + pybop.Observed("Current function [A]", solution["Current [A]"].data), + pybop.Observed("Voltage [V]", solution["Terminal voltage [V]"].data), ] # Define model From 4735367adcab4830dffca70eecaf1b48218e7522 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 18 Aug 2023 19:32:00 +0100 Subject: [PATCH 18/18] Updt. to default parameter set for a stable & working state --- tests/unit/test_parameterisation.py | 33 +++++++++++++++-------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/tests/unit/test_parameterisation.py b/tests/unit/test_parameterisation.py index a55dae54b..6bf4f6e92 100644 --- a/tests/unit/test_parameterisation.py +++ b/tests/unit/test_parameterisation.py @@ -1,38 +1,38 @@ import pybop -import pytest import pybamm +import pytest import numpy as np -import pandas as pd class TestParameterisation: - def getdata(self, x1, x2): + def getdata(self, x0): model = pybamm.lithium_ion.SPM() - params = pybamm.ParameterValues("Chen2020") + params = model.default_parameter_values params.update( { - "Negative electrode active material volume fraction": x1, - "Positive electrode active material volume fraction": x2, + "Negative electrode active material volume fraction": x0[0], + "Positive electrode active material volume fraction": x0[1], } ) experiment = pybamm.Experiment( [ ( - "Discharge at 0.5C for 5 minutes (1 second period)", + "Discharge at 2C for 5 minutes (1 second period)", "Rest for 2 minutes (1 second period)", - "Charge at 0.5C for 2.5 minutes (1 second period)", + "Charge at 1C for 5 minutes (1 second period)", "Rest for 2 minutes (1 second period)", ), ] * 2 ) sim = pybamm.Simulation(model, experiment=experiment, parameter_values=params) - return sim.solve(initial_soc=0.5) + return sim.solve() def test_rmse(self): # Form observations - solution = self.getdata(0.6, 0.6) + x0 = np.array([0.55, 0.63]) + solution = self.getdata(x0) observations = [ pybop.Observed("Time [s]", solution["Time [s]"].data), @@ -42,19 +42,19 @@ def test_rmse(self): # Define model model = pybop.models.lithium_ion.SPM() - model.parameter_set = pybop.ParameterSet("pybamm", "Chen2020") + model.parameter_set = model.pybamm_model.default_parameter_values # Fitting parameters params = [ pybop.Parameter( "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.1), - bounds=[0.05, 0.95], + prior=pybop.Gaussian(0.5, 0.05), + bounds=[0.35, 0.75], ), pybop.Parameter( "Positive electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.1), - bounds=[0.05, 0.95], + prior=pybop.Gaussian(0.65, 0.05), + bounds=[0.45, 0.85], ), ] @@ -67,4 +67,5 @@ def test_rmse(self): signal="Voltage [V]", method="nlopt" ) # Assertions - np.testing.assert_allclose(last_optim, 1e-1, atol=5e-2) + np.testing.assert_allclose(last_optim, 1e-3, atol=1e-2) + np.testing.assert_almost_equal(results, x0, decimal=1)