Plot capacity fade over throughput capacity

[Sunnyturicum]

Minimal working example

model = pybamm.lithium_ion.SPMe(
    {
        "lithium plating": "partially reversible",
        "lithium plating porosity change": "true",  # alias for "SEI porosity change"
    }
)

prada2013_params = pybamm.ParameterValues("OKane2022")

pybamm.set_logging_level("NOTICE")
cycle_number = 50
exp = pybamm.Experiment(
    [("Discharge at 1C until 2.5 V",  # ageing cycles
    "Charge at 0.3C until 3.65 V",
    "Hold at 3.65 V until C/20",
    "Rest for 4 hours",)] * cycle_number
)

temperatures = [273.15, 288.15] #-5, 5, 15°C, 14.3 268.15, 278.15, 287.45, 288.15

solutions = {} # dict, key is temp and values is solution object

for temp in temperatures:
    param = prada2013_params
    param.update({"Ambient temperature [K]": temp})

    sim = pybamm.Simulation(model, parameter_values=param, experiment=exp)
    safe_solver = pybamm.CasadiSolver(mode="safe", dt_max=120)
    sol = sim.solve(solver=safe_solver, calc_esoh=False)
    solutions[temp] = sol
    print(solutions[temp].data)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define the trend function with more descriptive variable names
def exponential_decay(x, amplitude, decay_rate, offset):
    # The trend function takes the inputs x, amplitude, decay_rate, and offset
    # and returns the value amplitude * exp(-decay_rate * x) + offset
    return amplitude * np.exp(-decay_rate * x) + offset

# Create a figure and axis for the plot
fig, ax = plt.subplots(figsize=(8, 6))
plt.rcParams.update({'font.size': 12})
line_colors = ["lightblue", "coral", "red", "cyan", "magenta", "yellow"]

# Loop through the solutions and their corresponding temperatures
for i, (temperature, sol) in enumerate(solutions.items()):
    throughput_capacity = sol["Throughput capacity [A.h]"].entries
    capacity_loss = sol["Loss of capacity to lithium plating [A.h]"].entries

    color = line_colors[i % len(line_colors)]
    ax.plot(
        throughput_capacity,
        capacity_loss,
        label="This work: " + str(temperature - 273.15) + " °C",
        color=color,
    )

    # Find the index of the minimum value in capacity_loss
    min_index = np.argmin(capacity_loss)
    
    # Fit the trendline through the minimum point using curve_fit
    initial_guess = [1, 1, min(capacity_loss)]  # Initial guess for parameters (amplitude, decay_rate, offset)
    fitted_params, _ = curve_fit(exponential_decay, throughput_capacity, capacity_loss, p0=initial_guess)
    
    # Generate x values for the trendline
    trendline_x = np.linspace(min(throughput_capacity), max(throughput_capacity), 100)
    
    # Calculate y values using the fitted parameters
    trendline_y = exponential_decay(trendline_x, *fitted_params)

    # Plot the trendline
    ax.plot(trendline_x, trendline_y, linestyle='dashed', color=color, alpha=0.7)

# Add labels and legend
ax.set_xlabel("Throughput capacity [A.h]", fontsize=14)
ax.set_ylabel("Capacity loss [A.h]", fontsize=14)
ax.legend()

# Adjust layout and display the plot
plt.tight_layout()
plt.show()

Problem

I want to plot the irreversible capacity loss that is steadily increasing at the bottom of the oscillating line. However, if I use the irreversible lithium plating model I overshoot. My approach was to plot a trendline. However, I think there must be a better way within Pybamm to plot this line.

[brosaplanella]

Have you tried using the summary variables (see this notebook)?

[Sunnyturicum]

Yes this would result into the following graph:

[rtimms]

You can get the solution broken down by cycles (solution.cycles) and steps (e.g. solution.cycles[0].steps) and then access variables from those. I think you are just wanting to plot the capacity loss at the end of a given step on each cycle, yes?

[ganmad]

Hi @Sunnyturicum
Have you managed to plot Capacity or Capacity loss vs Throughput capacity? Thanks.

[Sunnyturicum]

Hi @ganmad

This is what I did:

for cycle in cycle_numbers-1:
        throughput_cycle = sol.cycles[cycle]["Throughput capacity [A.h]"].data
        capacity_cycle = sol.cycles[cycle]["Loss of capacity to lithium plating [A.h]"].data
        min_index = np.argmin(capacity_cycle)

        capacity_value = capacity_cycle[min_index]
        corresponding_throughput_value = throughput_capacity[min_index]
        x_value.append(cycle_start_throughput + corresponding_throughput_value)
        y_value.append(capacity_value)
        cycle_start_throughput += throughput_cycle[-1]

Hopefully this helps.

[ganmad]

Thank you @Sunnyturicum. I will give it a try with your suggestion. Regards.