Designing degradation protocols¶
This notebook is the protocol catalogue: everything the two builders can express, and how to check a protocol before spending API time on it.
In this notebook
- The ageing protocols: voltage window, partial SoC window, calendar storage, and fully custom cycles
- Fast charge against the anode potential, worked end to end against a step charge
- Measured drive cycles imported from CSV
- The RPT in depth: capacity checks, HPPC pulse maps, rests, and an optional differential-analysis cycle.
import pandas as pd
import plotly.express as px
from breathe_simulate import api_interface as api
from breathe_simulate.ageing import AgeingCycler, RptCycler
cell_name = "Molicel P45B"
CAP_AH = 4.5
ageing_builder = AgeingCycler(selected_unit="C", cell_capacity=CAP_AH)
rpt_builder = RptCycler(selected_unit="C", cell_capacity=CAP_AH)
voltage_window = ageing_builder.cyclic(
I_chg=1.0,
I_dch=-1.0,
I_cut=0.05,
V_max=4.2,
V_min=2.5,
t_rest_s=300, # rest after each leg
t_max_cv_s=3600, # CV timeout
)
ageing_builder.plot_preview(voltage_window)
Previews are coulomb-count estimates, not simulations. Pass
soc_ref_capacity_ah=... to preview against a faded capacity.
1b. Partial SoC window¶
SoC bounds make the legs amp-hour limited, counted against the campaign's
SoC reference capacity. soc_reference="latest_rpt" (the default) re-bases
the window at each RPT, "bol" pins it to beginning of life.
A windowed leg ends at the first condition met: the coulomb-counted bound, the voltage limit, or a cutoff on the model's tracked SoC.
soc_window = ageing_builder.cyclic(
I_chg=0.33,
I_dch=-1.0,
I_cut=0.05,
V_max=4.2,
V_min=2.5,
soc_min=0.2,
soc_max=0.8,
)
# a faded cell shrinks the coulomb-counted window, so preview it both ways
ageing_builder.plot_preview(soc_window, soc_ref_capacity_ah=4.0)
1c. Calendar (storage) ageing¶
One storage block runs between RPTs. stop_on_time_days is the natural end
criterion. To alternate cycling and storage, pass a list of phases as the
ageing cycler.
calendar = ageing_builder.calendar(
storage_soc=0.5,
storage_days=30,
I_dch=-1.0,
I_chg=0.33,
I_cut=0.05,
V_max=4.2,
V_min=2.5,
)
ageing_builder.plot_preview(calendar)
1d. Fully custom cycles¶
custom() accepts what the run_sim CUSTOM cycler accepts: command
strings, structured step dictionaries, and drive-cycle arrays. Custom steps
carry their own units. Sections 2 and 3 use it.
2. Fast charge against the anode potential¶
Plating is governed by the anode potential versus Li/Li+ and becomes favourable below 0 mV. Measuring it needs a reference electrode, so no product cell reports it. The model computes it, so a charge can be controlled against it directly.
The alternative is a step charge: a current schedule that steps down as the cell reaches higher SoC.
step_charge = ageing_builder.custom(
experiment_text=[
"Charge at 13.5 A until 3.9 V", # 3C
"Charge at 9.0 A until 4.05 V", # 2C
"Charge at 4.5 A until 4.2 V", # 1C
"Hold at 4.2 V until 0.225 A",
"Rest for 600 s",
"Discharge at 4.5 A until 2.5 V",
"Rest for 600 s",
],
period="10 seconds",
)
The anode-potential version needs no schedule. Set the step charge's highest current as a ceiling and give the anode potential a floor. The current derates whenever holding the ceiling would cross the threshold.
ANODE_FLOOR_MV = 60.0 # the plating margin to hold
anode_limited = ageing_builder.custom(
experiment_text=[
{
"type": "vne_charge",
"current_a": 13.5, # 3C ceiling, derated as needed
"until_voltage_v": 4.2,
"anode_potential_threshold_mV": ANODE_FLOOR_MV,
},
"Hold at 4.2 V until 0.225 A",
"Rest for 600 s",
"Discharge at 4.5 A until 2.5 V",
"Rest for 600 s",
],
period="10 seconds",
)
# on the same step: temperature_threshold_degC derates to hold a temperature
# ceiling, and cut_off_current_a ends the step once the derated current
# falls below it
Now age both with identical campaign settings.
The headline metric is the time a 10 % to 80 % charge takes, over life. SoC
uses the test-bench convention: 100 % at the check-up's measured full,
counted against the latest measured capacity (soc_reference,
"latest_rpt" by default). A degraded cell has fewer amp hours between 10 %
and 80 %.
fast_charge_rpt = rpt_builder.build(
I_chg=1.0,
I_cut=0.05,
V_max=4.2,
V_min=2.5,
capacity_checks=[{"current": 1.0, "reference": True}],
pulses=[{"soc": 0.5, "current": -2.0, "duration_s": 30, "reference": True}],
t_rest_s=600,
t_equilibration_s=1800,
t_max_cv_s=3600,
reset_temperature=True,
)
CAMPAIGN = dict(
rpt_cycler=fast_charge_rpt,
rpt_every_n_cycles=50,
max_cycles=200,
initialTemperature_degC=25.0,
ambientTemperature_degC=25.0,
heatTransferCoefficient=35.0,
charge_time_soc_window=(0.1, 0.8),
return_charge_analysis=True,
)
stepped = api.run_ageing_sim(cell_name, step_charge, **CAMPAIGN)
guarded = api.run_ageing_sim(cell_name, anode_limited, **CAMPAIGN)
PROTOCOLS = {"step charge": stepped, f"anode floor {ANODE_FLOOR_MV:g} mV": guarded}
Running ageing campaign on 'Molicel P45B' (up to 200 cycles)... Campaign finished: max_cycles at cycle 200 (wall clock 00:01:31) Running ageing campaign on 'Molicel P45B' (up to 200 cycles)... Campaign finished: max_cycles at cycle 200 (wall clock 00:01:51)
2a. Speed and plating margin over life¶
. The charge-time column is named after its window (charge_time_column).
trend = pd.concat(
{name: r.cycles for name, r in PROTOCOLS.items()}, names=["protocol"]
).reset_index("protocol")
charge_col = stepped.charge_time_column # e.g. charge_time_10_80_pct_s
px.line(
trend,
x="cycle_number",
y=charge_col,
color="protocol",
title="Time to charge 10 % to 80 %, over life",
labels={charge_col: "Charge time [s]"},
)
px.line(
trend,
x="cycle_number",
y="min_anode_potential_mV",
color="protocol",
title="Minimum anode potential per cycle (0 mV = plating threshold)",
)
trend.groupby("protocol").agg(
charge_first_min=(charge_col, lambda s: s.iloc[0] / 60.0),
charge_last_min=(charge_col, lambda s: s.iloc[-1] / 60.0),
min_anode_potential_mV=("min_anode_potential_mV", "min"),
plating_lli_pct=("lli_plating [%]", "last"),
).round(3)
| charge_first_min | charge_last_min | min_anode_potential_mV | plating_lli_pct | |
|---|---|---|---|---|
| protocol | ||||
| anode floor 60 mV | 18.303 | 20.496 | 60.266 | 5.938 |
| step charge | 18.976 | 19.067 | 36.272 | 5.961 |
Fresh, the guarded charge is faster: it holds the 3C ceiling for as long as the anode allows, while the step charge drops to 2C at a preset voltage. It also holds a bigger plating margin.
Aged, they diverge. The step charge's time stays roughly flat because two effects cancel: the window is fewer amp hours, but higher resistance reaches each voltage stage earlier. Its margin erodes cycle by cycle. The guarded charge keeps the margin and slows down instead.
2b. Capacity, resistance and energy¶
Capacity and SoH come from the reference check, dcir_ohm from the
reference pulse, and the throughput columns accumulate over the campaign.
checkups = pd.concat({name: r.rpt for name, r in PROTOCOLS.items()}, names=["protocol"])
checkups[
[
"cycle_number",
"capacity_Ah",
"soh_pct",
"dcir_ohm",
"abs_ah_throughput_Ah",
"wh_in",
"wh_out",
"elapsed_time_s",
]
].round(4)
| cycle_number | capacity_Ah | soh_pct | dcir_ohm | abs_ah_throughput_Ah | wh_in | wh_out | elapsed_time_s | ||
|---|---|---|---|---|---|---|---|---|---|
| protocol | rpt_number | ||||||||
| step charge | 0 | 0 | 4.4376 | 100.0000 | 0.0145 | 19.9659 | 33.7155 | 39.2210 | 2.175084e+04 |
| 1 | 50 | 4.3330 | 97.6444 | 0.0153 | 475.6947 | 913.0228 | 858.0251 | 3.995779e+05 | |
| 2 | 100 | 4.2397 | 95.5406 | 0.0158 | 921.2990 | 1774.8830 | 1659.3059 | 7.725019e+05 | |
| 3 | 150 | 4.1512 | 93.5473 | 0.0161 | 1357.4731 | 2620.1743 | 2444.3285 | 1.140595e+06 | |
| 4 | 200 | 4.0663 | 91.6328 | 0.0164 | 1784.6446 | 3449.4863 | 3213.8307 | 1.504028e+06 | |
| anode floor 60 mV | 0 | 0 | 4.4376 | 100.0000 | 0.0145 | 19.9659 | 33.7155 | 39.2210 | 2.175084e+04 |
| 1 | 50 | 4.3332 | 97.6485 | 0.0153 | 475.0604 | 911.3262 | 858.0622 | 3.984743e+05 | |
| 2 | 100 | 4.2404 | 95.5576 | 0.0158 | 920.0694 | 1770.8680 | 1659.4319 | 7.755246e+05 | |
| 3 | 150 | 4.1526 | 93.5774 | 0.0161 | 1355.7112 | 2613.6355 | 2444.6289 | 1.150720e+06 | |
| 4 | 200 | 4.0682 | 91.6763 | 0.0164 | 1782.4077 | 3440.3471 | 3214.3938 | 1.523236e+06 |
# end of campaign, side by side
final = checkups.groupby("protocol").last()
pd.DataFrame(
{
"capacity_Ah": final["capacity_Ah"],
"SoH [%]": final["soh_pct"],
"DCIR [mOhm]": 1000 * final["dcir_ohm"],
"throughput [Ah]": final["abs_ah_throughput_Ah"],
"energy delivered [kWh]": final["wh_out"] / 1000.0,
"round-trip efficiency [%]": 100 * final["wh_out"] / final["wh_in"],
"duration [days]": final["elapsed_time_s"] / 86400.0,
}
).round(3)
| capacity_Ah | SoH [%] | DCIR [mOhm] | throughput [Ah] | energy delivered [kWh] | round-trip efficiency [%] | duration [days] | |
|---|---|---|---|---|---|---|---|
| protocol | |||||||
| anode floor 60 mV | 4.068 | 91.676 | 16.418 | 1782.408 | 3.214 | 93.432 | 17.630 |
| step charge | 4.066 | 91.633 | 16.402 | 1784.645 | 3.214 | 93.168 | 17.408 |
Both cells end up in a similar place: same voltage window, same discharge, so throughput, energy, fade and resistance growth all land close together. At these rates the guarded protocol buys its margin without a capacity or resistance penalty.
Equal cycle count is not equal calendar time.
duration [days] shows the difference. Normalising fade by energy delivered
and by time removes it.
fade = checkups.groupby("protocol").agg(
soh_lost_pct=("soh_pct", lambda s: s.iloc[0] - s.iloc[-1]),
kwh_delivered=("wh_out", lambda s: s.iloc[-1] / 1000.0),
days=("elapsed_time_s", lambda s: s.iloc[-1] / 86400.0),
)
fade["pct_per_kWh"] = fade["soh_lost_pct"] / fade["kwh_delivered"]
fade["pct_per_day"] = fade["soh_lost_pct"] / fade["days"]
fade.round(4)
| soh_lost_pct | kwh_delivered | days | pct_per_kWh | pct_per_day | |
|---|---|---|---|---|---|
| protocol | |||||
| anode floor 60 mV | 8.3237 | 3.2144 | 17.6300 | 2.5895 | 0.4721 |
| step charge | 8.3672 | 3.2138 | 17.4077 | 2.6035 | 0.4807 |
At 25 °C on this power cell the plating LLI barely differs between the two. The anode margin is comfortable at these rates. The floor matters where the margin is tight: cold charging, higher rates, or an aged cell.
2c. Where the charge-time numbers come from¶
return_charge_analysis=True shows the definition with the numbers, in
particular which SoC scale is used. The model's tracked SoC and the
campaign's protocol SoC differ once the cell has degraded.
guarded.charge_analysis_definition
{'soc_window': [0.1, 0.8],
'soc_scale': 'protocol (latest_rpt)',
'soc_source': "PROTOCOL SoC: 100 % at the RPT's defined full (the CC-CV anchor charge every check-up repeats identically), counted against the latest RPT's reference capacity, re-measured at every check-up, so the window's amp-hour span SHRINKS as the cell fades. At a fixed charging current the charge time therefore FALLS over life, exactly as it does on a test bench, until resistance growth (an earlier, longer CV phase) outweighs the shrinking window. Set soc_reference='bol' to pin the window to the baseline RPT's capacity instead: a fixed amp-hour span whose charge time grows with resistance.",
'soc_nominal_capacity_Ah': 4.4999999999999964,
'soc_protocol_reference_capacity_Ah': 4.068185839946252,
'soc_scale_caveat': "Every SoC value in this analysis is on the protocol scale above, the same scale the campaign's coulomb-counted set points (DCIR pulse placement, SoC-window ageing bounds) use. The model's internal tracked SoC keeps the nominal reference and is not what these columns report.",
'charge_time_definition': 'Wall-clock time between the interpolated instants at which the SoC rises through the low and then the high bound. The window is located as a RISE (down to the low bound, then up to the high bound), so the discharge that opens a conventional cycle is not mistaken for the charge and the rest at the bottom of the cycle is excluded. Anything inside the window counts, including a CV hold or a pause, because that is time a user waits at a charger.',
'not_traversed_behaviour': 'Cycles that never rise through the window (a narrower SoC window, calendar storage, a discharge-only cycle) report null rather than a guessed number.'}
The per-cycle table has the SoC at both boundaries, throughput, mean and peak current, boundary voltages, and the minimum anode potential and peak temperature inside the window. First row against last: same ceiling, lower mean C-rate, longer charge, margin held.
guarded.charge_analysis[
[
"cycle_number",
"charge_time_s",
"soc_at_window_entry",
"soc_at_window_exit",
"throughput_Ah",
"mean_c_rate",
"max_c_rate",
"min_anode_potential_mV",
"max_temperature_degC",
]
].iloc[[0, -1]].round(3)
| cycle_number | charge_time_s | soc_at_window_entry | soc_at_window_exit | throughput_Ah | mean_c_rate | max_c_rate | min_anode_potential_mV | max_temperature_degC | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1098.194 | 0.1 | 0.8 | 3.108 | 2.264 | 3.0 | 60.396 | 36.689 |
| 199 | 200 | 1229.750 | 0.1 | 0.8 | 2.909 | 1.892 | 3.0 | 60.288 | 35.935 |
3. Drive cycles from CSV¶
from_csv() turns a logged duty cycle into a repeatable ageing cycle. The
CSV needs a time column in seconds and a value column: current in amps
(discharge negative) or power in watts with type="power" (discharge
positive). With power input the solver delivers the commanded power at
whatever current the voltage requires. An optional ambient-temperature
column splits the profile into legs, each with its own ambient. pre steps
run before the profile every cycle.
Command strings accept watts too: "Discharge at 30 W until 2.5 V". Power
steps are validated against the calibration window using the worst-case
current, the power divided by the lower voltage cut-off. That is
conservative on purpose: the 3C discharge calibration allows 33.75 W at the
2.5 V floor, although 33.75 W at mid SoC only draws about 2C.
Give the ambient column as setpoints, not a raw sensor trace. Each contiguous setpoint becomes one simulated step, so a column that drifts sample by sample turns one drive cycle into hundreds of steps. A column that changes on almost every row is rejected.
In practice you would point from_csv at your own telemetry export. To keep
the example self-contained, the next cell synthesises a day of driving and
writes drive_cycle.csv: a cold urban start, a suburban link road, and a
hot highway leg with a 32 W merge and a regen exit.
import numpy as np
# one day of driving, synthesised as a stand-in for a real telemetry export.
# Power convention: watts, positive = discharge, negative = regen charging.
# Sized so the validator's worst case current (power / 2.5 V floor) stays
# inside this cell's 3C discharge and 3.4C charge calibration limits.
urban = []
for peak_w in (17.1, 23.5, 20.3, 27.7, 14.9, 25.6, 19.2, 21.3):
urban += [
*np.linspace(0.0, peak_w, 8), # pull away
*[peak_w] * 6, # accelerate
*np.linspace(peak_w, 3.2, 10), # back off to a coast
*[3.2] * 12, # coast
*np.linspace(3.2, -6.4, 5), # brake into regen
*[0.0] * 15, # stand at the lights
]
suburban = [
*np.linspace(0.0, 12.8, 12),
*[12.8] * 90, # link-road cruise
*np.linspace(12.8, -4.3, 8),
*[9.6] * 70, # dip, second cruise
*np.linspace(9.6, 0.0, 10),
]
highway = [
*np.linspace(0.0, 32.0, 10),
*[32.0] * 25, # hard merge
*(16.0 - 3.7 * np.sin(np.linspace(0, 12 * np.pi, 240))), # wavy cruise
*np.linspace(16.0, -16.0, 12),
*[0.0] * 20, # regen exit
]
power_w = np.round(np.concatenate([urban, suburban, highway]), 2)
ambient = [10.0] * len(urban) + [25.0] * len(suburban) + [35.0] * len(highway)
pd.DataFrame(
{
"time_s": np.arange(power_w.size, dtype=float),
"power_W": power_w,
"ambient_degC": ambient,
}
).to_csv("drive_cycle.csv", index=False)
drive = pd.read_csv("drive_cycle.csv")
print(f"{len(drive)} s, net {drive['power_W'].sum() / 3600:.2f} Wh drawn")
drive.head()
945 s, net 2.69 Wh drawn
| time_s | power_W | ambient_degC | |
|---|---|---|---|
| 0 | 0.0 | -0.00 | 10.0 |
| 1 | 1.0 | 2.44 | 10.0 |
| 2 | 2.0 | 4.88 | 10.0 |
| 3 | 3.0 | 7.31 | 10.0 |
| 4 | 4.0 | 9.75 | 10.0 |
px.line(
drive,
x="time_s",
y="power_W",
color="ambient_degC",
title="The measured duty cycle, coloured by ambient leg",
labels={"time_s": "Time [s]", "power_W": "Power [W]"},
)
ageing_drive = ageing_builder.from_csv(
"drive_cycle.csv",
time_col="time_s",
value_col="power_W",
type="power", # watts, discharge positive
ambient_col="ambient_degC",
# top up before each day, stopping short of full so the first regen leg
# has headroom instead of tripping the upper cut-off
pre=["Charge at 4.5 A until 4.1 V", "Rest for 600 s"],
post=["Rest for 600 s"],
# control-array legs carry no natural termination, so guard them: the leg
# ends gracefully at the first bound met instead of aborting the campaign
V_guard_min=2.6,
V_guard_max=4.15,
)
Power profiles cannot be previewed: a power leg's current depends on the
terminal voltage, which only a simulation knows, so plot_preview refuses.
The CSV plot above serves as the preview. Current profiles preview as usual.
Now simulate it. The timeseries reports the current the solver used to deliver each commanded watt, rising as the voltage sags.
drive_result = api.run_ageing_sim(
cell_name,
ageing_drive,
rpt_cycler=None,
max_cycles=5,
initialTemperature_degC=25.0,
ambientTemperature_degC=25.0,
heatTransferCoefficient=35.0,
return_timeseries=True,
timeseries_max_points_per_segment=120,
)
drive_result.plot_timeseries(["Current [A]", "Cell temperature [°C]"])
Running ageing campaign on 'Molicel P45B' (up to 5 cycles)... Campaign finished: max_cycles at cycle 5 (wall clock 00:00:29)
Every timeseries sample carries a step_number. For a drive cycle each item
is its own step: the pre steps, one step per ambient leg, then the post
steps. ageing_protocol lists what each index means.
drive_result.ageing_protocol
[{'phase': 0,
'steps': ['Charge at 4.5 A until 4.1 V',
'Rest for 600 s',
"drive-cycle leg: power array 'csv_profile_leg0' (ambient 10 degC)",
"drive-cycle leg: power array 'csv_profile_leg1' (ambient 25 degC)",
"drive-cycle leg: power array 'csv_profile_leg2' (ambient 35 degC)",
'Rest for 600 s']}]
Colouring one cycle's current by step_number separates the legs:
one_cycle = drive_result.timeseries.query(
"segment == 'ageing' and cycle_number == 5"
).assign(minutes=lambda d: (d["Time [s]"] - d["Time [s]"].min()) / 60.0)
px.line(
one_cycle,
x="minutes",
y="Current [A]",
color=one_cycle["step_number"].astype("Int64").astype(str),
title="Drive cycle current by step_number (ageing cycle 100)",
labels={"minutes": "Time into cycle [min]", "color": "step_number"},
)
4. The RPT in depth¶
The RPT measures the declared capacities, then the pulse resistances. All currents, voltages and rests are configurable. The step order is fixed so results stay comparable across campaigns and with lab RPTs. Each RPT starts from a conditioned, defined-empty state.
4a. Capacity checks and an HPPC pulse map¶
capacity_checks and pulses are lists, each with one reference entry.
Several capacity checks give a rate-dependent capacity map. Pulse SoC set
points are coulomb-counted against the reference capacity measured in the
same RPT, so "50 % SoC" tracks the degraded cell.
rpt_hppc = rpt_builder.build(
I_chg=1.0,
I_cut=0.05,
V_max=4.2,
V_min=2.5,
capacity_checks=[
{"current": 1.0, "reference": True}, # SoH and the SoC anchor
{"current": 0.33}, # slow check, for the rate gap
],
t_equilibration_s=1800,
pulses=[ # signed C-rates: negative is a discharge pulse
{"soc": 0.8, "current": -2.0, "duration_s": 30},
{"soc": 0.8, "current": 1.0, "duration_s": 30},
{"soc": 0.5, "current": -2.0, "duration_s": 30, "reference": True},
{"soc": 0.5, "current": 1.0, "duration_s": 30},
{"soc": 0.2, "current": -2.0, "duration_s": 30},
{"soc": 0.2, "current": 1.0, "duration_s": 30},
],
t_rest_s=300,
t_max_cv_s=3600,
reset_temperature=True, # run the RPT at the reference temperature
)
rpt_builder.plot_preview(rpt_hppc)