3. Cycling data

In this example, we show how to use full-cell cycling data to fit model parameters. Specifically, we will fit the negative electrode reference exchange-current density. The workflow shown in this example is very general and can take any cycling data and fit the model to any collected data. Typically we fit the voltage and/or temperature.

Before running this example, make sure to run the script true_parameters/generate_data.py to generate the synthetic data.

from pathlib import Path

import ionworkspipeline as iwp
import iwutil
import pybamm
from true_parameters.parameters import full_cell

Load the data

This workflow can fit across multiple data at once. Here we will fit the time vs voltage data collected from three constant current experiments. We load the data into a dictionary whose keys are the C-rates.

data = {}
for rate in [0.1, 0.5, 1]:
    data[f"{rate} C"] = iwutil.read_df(
        Path("synthetic_data") / "full_cell" / f"{str(rate).replace('.', '_')}C.csv"
    )

Set up the parameters to fit

Now we set up our initial guesses and bounds for the parameters. As in the previous example, we create a dictionary of Parameter objects. We assume a symmetric Butler-Volmer exchange-current density. We create a function for the exchange-current density and fit a scalar reference exchange-current density.

def j0_n(c_e, c_s_surf, c_s_max, T):
    j0_ref = pybamm.Parameter(
        "Negative electrode reference exchange-current density [A.m-2]"
    )
    alpha = 0.5
    c_e_ref = 1000

    return (
        j0_ref
        * (c_e / c_e_ref) ** (1 - alpha)
        * (c_s_surf / c_s_max) ** alpha
        * (1 - c_s_surf / c_s_max) ** (1 - alpha)
    )


parameters = {
    "Negative electrode exchange-current density [A.m-2]": j0_n,
    "Negative electrode reference exchange-current density [A.m-2]": iwp.Parameter(
        "j0_n [A.m-2]",
        initial_value=8,
        bounds=(0.1, 10),
    ),
}

Get the known parameters

In this example, we assume we already know all of the model parameters except for the negative electrode diffusivity and reference exchange-current density. We load the known parameters from the true parameters script.

true_params = full_cell()
known_params = {k: v for k, v in true_params.items() if k not in parameters}

Run the fit

Now we are ready to run the fit. We create the objective(s), the DataFit object, and run the fit.

We use the CurrentDriven objective function, which takes the input current from the data and fits the model output to the objective variables specified by the user. By default it fits the voltage. Since we are fitting multiple data sets, we create a corresponding dictionary of Objective objects, and the fit simultaneously optimizes across all data.

We pass a pybamm model to the objective function. Here we use the SPMe model, with the events removed to avoid triggering the voltage cut-off event. We also use the parameters keyword argument, which allows us to pass additional parameters to the objective function that are specific to each objective. Here we use it to match the initial voltage in the data.

Finally, we create the DataFit object and run the fit.

model = pybamm.lithium_ion.SPMe()
model.events = []
options = {"model": model}


def get_objective_parameters(data):
    V_initial = data["Voltage [V]"].iloc[0]
    return {
        "Initial voltage [V]": V_initial,
    }


# Create a dictionary of objectives, one for each data set
objectives = {}
for this_name, this_data in data.items():
    objectives[this_name] = iwp.objectives.CurrentDriven(
        this_data, options=options, parameters=get_objective_parameters(this_data)
    )

# Create the DataFit object and run the fit
cost = iwp.costs.RMSE()
optimizer = iwp.optimizers.ScipyMinimize(method="Nelder-Mead")
datafit = iwp.DataFit(objectives, parameters=parameters, cost=cost, optimizer=optimizer)
results = datafit.run(known_params)

We can plot the trace of the cost function and parameter values during the fit, and the fitted voltage curves.

_ = datafit.plot_trace()
_ = datafit.plot_fit_results()

We can take a look at the fitted parameters by accessing the parameter_values attribute of the results object.

results.parameter_values

{'Negative electrode exchange-current density [A.m-2]': <function __main__.j0_n(c_e, c_s_surf, c_s_max, T)>,
 'Negative electrode reference exchange-current density [A.m-2]': 0.6757812499999936}