Discrepancy in Degradation Calculations When Using Experiment Blocks

Hello everyone,
I have encountered an issue with simulating degradation mechanisms in PyBaMM when splitting a long experiment into multiple blocks. Here are the details:

I have real-world power data over a month, provided in 10-minute intervals. These data include charge and discharge cycles, which I have converted into experiment steps (experiment_steps). This results in:

  • Per day: 144 steps
  • Per month: 4320 steps

When attempting to run such a long experiment in a single simulation, I encounter an error, which I interpret as a timeout. To address this, I attempted to split the experiment into smaller blocks. After each block, I calculated the SEI thickness and nominal cell capacity at the end of the block and passed them as updated parameters to the next block.

To verify the validity of this approach, I created an experiment with 144 steps and executed it in two ways:

  1. Without splitting into blocks – a single experiment with all 144 steps.
  2. With blocks – dividing the 144 steps into 24 blocks of 6 steps each.

Observation:
When comparing the results, I noticed that the total degradation in A.h. was significantly higher in the block-wise execution compared to the single-run experiment.

Questions:

  1. Which variables and parameter values need to be correctly updated after each block to ensure the results are consistent with running the experiment as a single long simulation?
  2. Are there known limitations or challenges when simulating long experiments in PyBaMM that are split into smaller blocks?
  3. Are there best practices for handling such long experiments, particularly regarding numerical consistency when passing state variables (e.g., SEI thickness, capacity) between blocks?

Note:
The relevant code is attached to this discussion for reference.

Thank you in advance for any help or suggestions!


import pybamm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from battery_classes import Cell, Battery
import time
import json
import os

os.chdir(os.path.dirname(__file__))

# Load JSON
with open("experiment_steps.json", "r") as json_file:
    experiment_steps = json.load(json_file)

battery = Battery()
soc = battery.initial_soc

# Initialize
parameter_values = pybamm.ParameterValues("OKane2022")
model = pybamm.lithium_ion.SPM(
    {
        "SEI": "solvent-diffusion limited",
        "SEI porosity change": "false",
        "lithium plating": "partially reversible",
        "lithium plating porosity change": "false",  # alias for "SEI porosity change"
        "particle mechanics": ("swelling and cracking", "swelling only"),
        "SEI on cracks": "true",
        "loss of active material": "stress-driven",
        "calculate discharge energy": "true",  # for compatibility with older PyBaMM versions
    }
)
solver = pybamm.CasadiSolver(mode="fast")

# Initialize the simulation
block_size = 6  # Number of steps per block
test_steps = experiment_steps[:144]  # 144 steps = 1 day
total_capacity = 5.0

# Split the experiment into blocks
experiment_blocks = [test_steps[i:i + block_size] for i in range(0, len(test_steps), block_size)]

for block_index, block in enumerate(experiment_blocks):
    print(f"Solving block {block_index + 1} of {len(experiment_blocks)}...")
    print(block)
    experiment_deg = pybamm.Experiment(block, period="10 minutes")
    sim = pybamm.Simulation(model, experiment=experiment_deg, parameter_values=parameter_values, solver=solver)
    solution = sim.solve(initial_soc=soc)  

    # Evaluate results
    Q_max = parameter_values["Nominal cell capacity [A.h]"]  # Maximum capacity
    Q_current = solution["Discharge capacity [A.h]"].entries[-1]  # Discharge capacity
    soc_diff = Q_current / Q_max  # Recalculate SOC
    soc -= soc_diff   

    Q_SEI = solution["Loss of capacity to negative SEI [A.h]"].entries[-1]
    Q_SEI_cr = solution["Loss of capacity to negative SEI on cracks [A.h]"].entries[-1]
    Q_plating = solution["Loss of capacity to negative lithium plating [A.h]"].entries[-1]
    Q_side = solution["Total capacity lost to side reactions [A.h]"].entries[-1]
    Q_LLI = (solution["Total lithium lost [mol]"].entries[-1] * 96485.3 / 3600)  # Convert from mol to A.h
    print("Loss of capacity to negative SEI [A.h]", Q_SEI)
    print("Loss of capacity to negative SEI on cracks [A.h]", Q_SEI_cr)
    print("Loss of capacity to negative lithium plating [A.h]", Q_plating)
    print("Total capacity lost to side reactions [A.h]", Q_side)
    print("Total lithium lost [mol]", Q_LLI)

    # Calculate total capacity losses
    total_capacity_loss = Q_SEI + Q_SEI_cr + Q_plating + Q_side + Q_LLI    
    total_capacity -= total_capacity_loss
    print(f"Updated total capacity: {total_capacity:.4f} A.h")

    sei_avg_neg_in = solution["X-averaged negative inner SEI thickness [m]"].entries[-1]
    sei_avg_neg_out = solution["X-averaged negative outer SEI thickness [m]"].entries[-1]
    parameter_values.update({"Nominal cell capacity [A.h]": total_capacity, 
                             "Initial inner SEI thickness [m]": sei_avg_neg_in,
                             "Initial outer SEI thickness [m]": sei_avg_neg_out,
                             })

    print(f"Block {block_index + 1} successfully completed.")