Combining Objectives example

This notebook shows how to construct and combine custom Objective classes to fit multiple objectives at once (e.g. simultaneously fit two data sets). For design principles, see the user_guide. For detailed documentation on the Objective class, see the API reference.

import pybamm
import ionworkspipeline as iwp
import numpy as np
import pandas as pd
import pybammeis

First, we generate timeseries data for a CurrentDriven objective using the SPM with the differential surface form option enabled.

# Load the model and parameter values
model = pybamm.lithium_ion.SPM(options={"surface form": "differential"})
parameter_values = iwp.ParameterValues("Chen2020")

parameter_values.update(
    {
        "Negative electrode active material volume fraction": 0.75,
        "Positive electrode active material volume fraction": 0.665,
    }
)


def j0(c_e, c_s_surf, c_s_max, T):
    j0_ref = pybamm.Parameter(
        "Positive electrode reference exchange-current density [A.m-2]"
    )
    c_e_init = pybamm.Parameter("Initial concentration in electrolyte [mol.m-3]")

    return (
        j0_ref
        * (c_e / c_e_init) ** 0.5
        * (c_s_surf / c_s_max) ** 0.5
        * (1 - c_s_surf / c_s_max) ** 0.5
    )


parameter_values.update(
    {
        "Positive electrode reference exchange-current density [A.m-2]": 2,
        "Positive electrode exchange-current density [A.m-2]": j0,
        "Positive electrode double-layer capacity [F.m-2]": 0.3,
    },
    check_already_exists=False,
)

# Generate the data
t = np.arange(0, 3000, 3)
sim = iwp.Simulation(
    model, parameter_values=parameter_values, t_eval=[t[0], t[-1]], t_interp=t
)
sol = sim.solve()
data = pd.DataFrame(
    {x: sim.solution[x](t) for x in ["Time [s]", "Current [A]", "Voltage [V]"]}
)
# add noise to the data
sigma = 0.001

# Set seed for reproducibility
np.random.seed(0)
data["Voltage [V]"] += np.random.normal(0, sigma, len(data))

objective_current_driven = iwp.objectives.CurrentDriven(data, options={"model": model})

Next, we use the same model and parameter_values to create synthetic EIS data and the EIS objective.

# Create simulation
eis_sim = pybammeis.EISSimulation(model, parameter_values=parameter_values)

# Choose frequencies and calculate impedance
frequencies = np.logspace(-4, 4, 30)
sol = eis_sim.solve(frequencies)

# Store data as arrays of real and complex impedance
Z_Re = sol.real
Z_Im = sol.imag
data = pd.DataFrame(
    {"Frequency [Hz]": frequencies, "Z_Re [Ohm]": Z_Re, "Z_Im [Ohm]": Z_Im}
)

objective_eis = iwp.objectives.EIS(data, options={"model": model})

# Generate a Nyquist plot
_ = eis_sim.nyquist_plot()

Now, we define the fitting parameters that will be optimized using both the CurrentDriven and EIS data. We perform a Log10 transform on the cathode exchange-current density and double-layer capacity because these parameters span several orders of magnitude.

parameters = {
    "Negative electrode active material volume fraction": iwp.Parameter(
        "Anode active material frac.",
        bounds=(0.3, 0.8),
    ),
    "Positive electrode active material volume fraction": iwp.Parameter(
        "Cathode active material frac.",
        bounds=(0.3, 0.8),
    ),
    "Positive electrode reference exchange-current density [A.m-2]": iwp.transforms.Log10(
        iwp.Parameter(
            "j0_p_ref",
            bounds=(0.1, 10),
        )
    ),
    "Positive electrode double-layer capacity [F.m-2]": iwp.transforms.Log10(
        iwp.Parameter(
            "C_dl_p",
            bounds=(0.01, 1),
        )
    ),
}

Now, we combine the two objectives to optimize them simultaneously.

objectives = {
    "CurrentDriven": objective_current_driven,
    "EIS": objective_eis,
}

All Objectives have a default optimizer associated with them that is tuned for the particular usecase. A specific optimizer can always be given to the DataFit class, but if none is provided, the optimizer will default to the one given by the first objective.

The default optimizers for the EIS and CurrentDriven objectives are:

for name, objective in objectives.items():
    print(f"{name}: {objective.default_optimizer.name}")

CurrentDriven: SciPy Minimize Optimizer
EIS: SciPy Least Squares Optimizer

Next, we create the DataFit class combining both objectives. The DataFit uses the SciPy Minimize Optimizer provided by objective_current_driven.

data_fit = iwp.DataFit(objectives, parameters=parameters)

print(f"DataFit optimizer: {data_fit.optimizer.name}")

DataFit optimizer: SciPy Minimize Optimizer

Next, we run the data fit, passing the parameters that are not being fit as a dictionary.

params_for_pipeline = {k: v for k, v in parameter_values.items() if k not in parameters}
results = data_fit.run(params_for_pipeline)
for k, v in results.items():
    print(f"{k}: {parameter_values[k]:.3e} (true) {v:.3e} (fit)")

Negative electrode active material volume fraction: 7.500e-01 (true) 7.502e-01 (fit)
Positive electrode active material volume fraction: 6.650e-01 (true) 6.649e-01 (fit)
Positive electrode reference exchange-current density [A.m-2]: 2.000e+00 (true) 1.999e+00 (fit)
Positive electrode double-layer capacity [F.m-2]: 3.000e-01 (true) 2.966e-01 (fit)

Finally, we plot the results and the trace of the cost function and parameter values for both the CurrentDriven and EIS objectives.

_ = data_fit.plot_fit_results()
fig, axes = data_fit.plot_trace()

sorted_keys = np.asarray(list(parameters.keys()))[
    np.argsort([p.name for p in parameters.values()])
]

for ax, k in zip(axes[1:], sorted_keys):
    v = parameter_values[k]
    if data_fit.parameters_original[k].is_transformed:
        v = data_fit.parameters_original[k].transform(v)
    ax.axhline(v, color="red", linestyle="--", label="True")
    ax.legend()