# SOH with Drive Cycle Profiles with Coupled Degradation Models

**URL:** <https://pybamm.discourse.group/t/soh-with-drive-cycle-profiles-with-coupled-degradation-models/383>\
**Category:** PyBaMM Q&A\
**Created:** [18 August 2025 16:19 UTC](https://pybamm.discourse.group/t/soh-with-drive-cycle-profiles-with-coupled-degradation-models/383 "2025-08-18T16:19:36Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![AkilaJ](https://yyz2.discourse-cdn.com/free1/user_avatar/pybamm.discourse.group/akilaj/32/12_2.png) [@AkilaJ](https://pybamm.discourse.group/u/AkilaJ)\
**Post date:** [18 August 2025 16:19 UTC](https://pybamm.discourse.group/t/soh-with-drive-cycle-profiles-with-coupled-degradation-models/383/1 "2025-08-18T16:19:36Z")

</div>

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;

1. Why is the capacity too small
2. After the first cycle, why did the capacity suddenly drop?

Your suggestion will be really helpful, and thank you for your consideration.

 ![image](https://global.discourse-cdn.com/free1/uploads/pybamm/original/1X/030733f0a5bf7d9c53bf45b9093d165a264a0bca.png)

* * *

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()