Hi All and @DrSOKane
I am trying to use the coupled degradation model for drive cycle simulations. As a start, I wanted to use the US06 drive cycle with the coupled degradation models and obtain capacity (or SOH) variation against cycle number. However, I have obtained the abnormal behaviour, as can be seen below.
I have two questions;
- Why is the capacity too small
- After the first cycle, why did the capacity suddenly drop?
Your suggestion will be really helpful, and thank you for your consideration.
Code used to generate this plot.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pybamm
from pathlib import Path
-------------------------------
1) Model + parameters
-------------------------------
model = pybamm.lithium_ion.SPMe(
{
“SEI”: “solvent-diffusion limited”,
“SEI porosity change”: “true”,
“lithium plating”: “partially reversible”,
“lithium plating porosity change”: “true”,
“particle mechanics”: “swelling and cracking”,
“SEI on cracks”: “true”,
“loss of active material”: “stress-driven”,
“thermal”: “lumped”, # optional
}
)
param = pybamm.ParameterValues(“OKane2022”)
-------------------------------
2) Load US06 (time[s], current[A])
-------------------------------
CSV_PATH = Path(“US06.csv”)
if CSV_PATH.exists():
base_dc = pd.read_csv(CSV_PATH, comment=“#”, header=None).to_numpy()
else:
data_loader = pybamm.DataLoader()
base_dc = pd.read_csv(
data_loader.get_data(“US06.csv”), comment=“#”, header=None
).to_numpy()
assert base_dc.shape[1] == 2, “CSV must have two columns: time(s), current(A)”
base_dc[:, 0] -= base_dc[0, 0] # normalize start time
Safety: drop any duplicate timestamps in the base cycle
_, uniq_idx = np.unique(base_dc[:, 0], return_index=True)
base_dc = base_dc[np.sort(uniq_idx)]
T_base = float(base_dc[-1, 0]) # base cycle duration (s)
-------------------------------
3) Repeat cycle WITHOUT duplicate boundary times
Keep the first sample of every cycle, drop the last sample of each
block except the final one.
-------------------------------
def repeat_drive_cycle_no_overlap(dc: np.ndarray, repeats: int = 10) → np.ndarray:
t, i = dc[:, 0], dc[:, 1]
T = t[-1]
blocks =
for k in range(repeats):
block = np.column_stack([t + k * T, i])
if k < repeats - 1:
block = block[:-1] # drop last row to avoid dup with next block’s first row
blocks.append(block)
out = np.vstack(blocks)
Final guard: ensure strictly increasing time
assert np.all(np.diff(out[:, 0]) > 0), “t_eval must be strictly increasing”
return out
repeats = 10 # change as you like
drive_cycle = repeat_drive_cycle_no_overlap(base_dc, repeats=repeats)
-------------------------------
4) Apply as current input & simulate
-------------------------------
current_interp = pybamm.Interpolant(drive_cycle[:, 0], drive_cycle[:, 1], pybamm.t)
param[“Current function [A]”] = current_interp
solver = pybamm.CasadiSolver(mode=“fast”)
sim = pybamm.Simulation(model, parameter_values=param, solver=solver)
t_eval = drive_cycle[:, 0] # strictly increasing timestamps
solution = sim.solve(t_eval=t_eval, t_interp=“linear”)
-------------------------------
5) Extract signals
-------------------------------
time = solution[“Time [s]”].data
current = solution[“Current [A]”].data
Qd = solution[“Discharge capacity [A.h]”].data # cumulative discharge capacity
-------------------------------
6) Incremental discharge capacity per US06 window
-------------------------------
discharge_capacities =
cycle_numbers =
Qd_prev = 0.0
for k in range(1, repeats + 1):
t_lo, t_hi = (k - 1) * T_base, k * T_base
idx = np.where((time > t_lo) & (time <= t_hi))[0]
if idx.size == 0:
discharge_capacities.append(0.0)
cycle_numbers.append(k)
continue
Qd_max_k = float(np.max(Qd[idx])) # cumulative max within this window
cap_k = Qd_max_k - Qd_prev # incremental capacity for this cycle
discharge_capacities.append(cap_k)
cycle_numbers.append(k)
Qd_prev = Qd_max_k
-------------------------------
7) SOH vs cycle (relative to cycle 1)
-------------------------------
initial_capacity = discharge_capacities[0] if discharge_capacities else 1.0
soh = [100.0 * cap / initial_capacity for cap in discharge_capacities]
-------------------------------
8) Plots
-------------------------------
plt.figure(figsize=(8, 4.5))
plt.plot(cycle_numbers, discharge_capacities, “bo-”, label=“Discharge Capacity per cycle”)
plt.xlabel(“Cycle Number”)
plt.ylabel(“Capacity (A·h)”)
plt.title(“Capacity vs Cycle (US06)”)
plt.grid(True)
plt.legend()
plt.figure(figsize=(8, 4.5))
plt.plot(cycle_numbers, soh, “rs-”, label=“SOH”)
plt.xlabel(“Cycle Number”)
plt.ylabel(“SOH (%)”)
plt.title(“SOH vs Cycle (US06)”)
plt.grid(True)
plt.legend()
plt.show()