ECM generation¶

generate_ecm runs a virtual pulse characterisation on the physics model and fits an RC network to it, so the equivalent circuit model inherits the physics of the exact design you are working on.

1. A fresh ECM for a design¶

ecm_options sets up the fit. A run_sim style dict passed as design=... characterises a modified cell instead of the library one.

In [1]:

from breathe_simulate import api_interface as api
from breathe_simulate.ecm import ecm_options

cell_name = "Molicel P45B"

options = ecm_options(
    n_rc=2,
    soc_grid=[0.1, 0.3, 0.5, 0.7, 0.9],
    temp_grid_degC=[10.0, 25.0, 40.0],
    charge_c_rates=[1.0, 2.0],
    discharge_c_rates=[1.0, 3.0],
)

fresh = api.generate_ecm(cell_name, ecm=options)
fresh

from breathe_simulate import api_interface as api from breathe_simulate.ecm import ecm_options cell_name = "Molicel P45B" options = ecm_options( n_rc=2, soc_grid=[0.1, 0.3, 0.5, 0.7, 0.9], temp_grid_degC=[10.0, 25.0, 40.0], charge_c_rates=[1.0, 2.0], discharge_c_rates=[1.0, 3.0], ) fresh = api.generate_ecm(cell_name, ecm=options) fresh

Out[1]:

EcmResults(2RC, 5 SoC x 3 T x 2 dch / 2 chg rates)

The parameter tables¶

to_dataframe() flattens the fitted tables into one long table.

In [2]:

table = fresh.to_dataframe()
table.head(12)

table = fresh.to_dataframe() table.head(12)

Out[2]:

	soc	temperature_degC	c_rate	direction	R0_Ohm	R1_Ohm	C1_F	R2_Ohm	C2_F
0	0.1	10.0	1.0	discharge	0.009971	0.010867	1629.678595	0.022025	7323.504658
1	0.1	10.0	3.0	discharge	0.009971	0.010867	1629.678595	0.022025	7323.504658
2	0.1	25.0	1.0	discharge	0.007971	0.006179	1724.726773	0.008665	11601.253334
3	0.1	25.0	3.0	discharge	0.008105	0.007536	1110.445497	0.008687	8973.488684
4	0.1	40.0	1.0	discharge	0.007070	0.004223	1565.663322	0.004121	17473.713909
5	0.1	40.0	3.0	discharge	0.007186	0.004597	1225.185020	0.003470	17824.631275
6	0.3	10.0	1.0	discharge	0.008181	0.006393	2261.469746	0.011410	17831.717000
7	0.3	10.0	3.0	discharge	0.008156	0.005407	3175.427032	0.010981	16092.795338
8	0.3	25.0	1.0	discharge	0.006488	0.003984	2331.950971	0.005942	22394.349336
9	0.3	25.0	3.0	discharge	0.006588	0.003719	2582.923262	0.004840	23896.444299
10	0.3	40.0	1.0	discharge	0.005753	0.002933	2132.951552	0.003816	24039.072258
11	0.3	40.0	3.0	discharge	0.005843	0.002776	2219.873191	0.002846	27787.915255

r0, r(k) and c(k) return the raw 3D arrays over (SoC x temperature x C rate).

In [3]:

print("R0 discharge table shape:", fresh.r0("discharge").shape)
print("SoC grid:", fresh.soc_grid)
print("Temperature grid [degC]:", fresh.temp_grid_degC)
print("Capacity [Ah]:", fresh.capacity)

fresh.plot_parameters("R0", direction="discharge", c_rate=1.0)

print("R0 discharge table shape:", fresh.r0("discharge").shape) print("SoC grid:", fresh.soc_grid) print("Temperature grid [degC]:", fresh.temp_grid_degC) print("Capacity [Ah]:", fresh.capacity) fresh.plot_parameters("R0", direction="discharge", c_rate=1.0)

R0 discharge table shape: (5, 3, 2)
SoC grid: [0.1 0.3 0.5 0.7 0.9]
Temperature grid [degC]: [10. 25. 40.]
Capacity [Ah]: 4.4999999999999964

Each direction and each C rate gets its own table. Comparing them shows the asymmetry between charge and discharge pulses and how the resistance shifts with current.

In [4]:

import plotly.graph_objects as go

temp_index = 1  # 25 degC in the grid above

fig = go.Figure()
for direction, rates in (
    ("charge", fresh.charge_c_rates),
    ("discharge", fresh.discharge_c_rates),
):
    for i, rate in enumerate(rates):
        fig.add_trace(
            go.Scatter(
                x=fresh.soc_grid,
                y=1000 * fresh.r0(direction)[:, temp_index, i],
                mode="lines+markers",
                name=f"{rate:g}C {direction}",
            )
        )
fig.update_layout(
    title="R0 vs SoC at 25 degC, both directions and pulse currents",
    xaxis_title="SoC",
    yaxis_title="R0 [mOhm]",
    legend_title="Pulse",
)
fig

import plotly.graph_objects as go temp_index = 1 # 25 degC in the grid above fig = go.Figure() for direction, rates in ( ("charge", fresh.charge_c_rates), ("discharge", fresh.discharge_c_rates), ): for i, rate in enumerate(rates): fig.add_trace( go.Scatter( x=fresh.soc_grid, y=1000 * fresh.r0(direction)[:, temp_index, i], mode="lines+markers", name=f"{rate:g}C {direction}", ) ) fig.update_layout( title="R0 vs SoC at 25 degC, both directions and pulse currents", xaxis_title="SoC", yaxis_title="R0 [mOhm]", legend_title="Pulse", ) fig

OCV, thermal constants and saving¶

The result also carries the OCV curves and thermal constants an ECM runtime needs, so it is a complete model on its own. The charge and discharge branches of the OCV show the hysteresis the mean curve averages out.

In [5]:

fresh.ocv.head()

fresh.ocv.head()

Out[5]:

	soc	mean_voltage_V	hysteresis_voltage_V	charge_voltage_V	discharge_voltage_V
0	0.000	2.663782	0.163802	2.827584	2.499981
1	0.001	2.685911	0.165880	2.851792	2.520031
2	0.002	2.705992	0.166806	2.872798	2.539186
3	0.003	2.724833	0.167096	2.891929	2.557738
4	0.004	2.742614	0.166951	2.909565	2.575664

In [6]:

fig = go.Figure()
for column, label in (
    ("charge_voltage_V", "charge branch"),
    ("discharge_voltage_V", "discharge branch"),
    ("mean_voltage_V", "mean"),
):
    fig.add_trace(
        go.Scatter(x=fresh.ocv["soc"], y=fresh.ocv[column], mode="lines", name=label)
    )
fig.update_layout(
    title="Open circuit voltage and its hysteresis",
    xaxis_title="SoC",
    yaxis_title="Voltage [V]",
    legend_title="Curve",
)
fig

fig = go.Figure() for column, label in ( ("charge_voltage_V", "charge branch"), ("discharge_voltage_V", "discharge branch"), ("mean_voltage_V", "mean"), ): fig.add_trace( go.Scatter(x=fresh.ocv["soc"], y=fresh.ocv[column], mode="lines", name=label) ) fig.update_layout( title="Open circuit voltage and its hysteresis", xaxis_title="SoC", yaxis_title="Voltage [V]", legend_title="Curve", ) fig

In [7]:

{
    key: f"curve with {len(value)} points" if isinstance(value, list) else value
    for key, value in fresh.thermal.items()
}

{ key: f"curve with {len(value)} points" if isinstance(value, list) else value for key, value in fresh.thermal.items() }

Out[7]:

{'cell_mass_kg': 0.06822167545418895,
 'cell_heat_capacity_J_kgK': 1200.0,
 'cell_surface_area_m2': 0.005370967488621659,
 'soc_dudt': 'curve with 1001 points',
 'dudt_VK': 'curve with 1001 points',
 'thermal_conductivity_in_w_mk': 42.23552620603363,
 'thermal_conductivity_through_w_mk': 1.0855399470076916}

In [8]:

from breathe_simulate.results import EcmResults

fresh.save("fresh_ecm.json")
reloaded = EcmResults.load("fresh_ecm.json")
reloaded

from breathe_simulate.results import EcmResults fresh.save("fresh_ecm.json") reloaded = EcmResults.load("fresh_ecm.json") reloaded

Out[8]:

EcmResults(2RC, 5 SoC x 3 T x 2 dch / 2 chg rates)

2. An aged ECM from a degradation campaign¶

Pass ecm=True (or an ecm_options(...) block) to run_ageing_sim and the same fit runs on the aged cell at the end of the campaign, with the final degradation applied to the physics. A lighter grid keeps the runtime down.

In [9]:

from breathe_simulate.ageing import AgeingCycler, RptCycler

CAP_AH = 4.5

ageing = AgeingCycler(selected_unit="C", cell_capacity=CAP_AH).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,
    t_max_cv_s=3600,
)
rpt = RptCycler(selected_unit="C", cell_capacity=CAP_AH).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": -1.0, "duration_s": 30, "reference": True}],
    t_equilibration_s=1800,
)

ageing_ecm = ecm_options(
    n_rc=2,
    soc_grid=[0.1, 0.5, 0.9],
    temp_grid_degC=[10.0, 25.0],
    charge_c_rates=[1.0],
    discharge_c_rates=[1.0, 3.0],
)

result = api.run_ageing_sim(
    cell_name,
    ageing,
    rpt_cycler=rpt,
    max_cycles=200,
    rpt_every_n_cycles=50,
    ecm=ageing_ecm,
)

aged = result.ecm
aged

from breathe_simulate.ageing import AgeingCycler, RptCycler CAP_AH = 4.5 ageing = AgeingCycler(selected_unit="C", cell_capacity=CAP_AH).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, t_max_cv_s=3600, ) rpt = RptCycler(selected_unit="C", cell_capacity=CAP_AH).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": -1.0, "duration_s": 30, "reference": True}], t_equilibration_s=1800, ) ageing_ecm = ecm_options( n_rc=2, soc_grid=[0.1, 0.5, 0.9], temp_grid_degC=[10.0, 25.0], charge_c_rates=[1.0], discharge_c_rates=[1.0, 3.0], ) result = api.run_ageing_sim( cell_name, ageing, rpt_cycler=rpt, max_cycles=200, rpt_every_n_cycles=50, ecm=ageing_ecm, ) aged = result.ecm aged

Running ageing campaign on 'Molicel P45B' (up to 200 cycles)...
Campaign finished: max_cycles at cycle 200 (wall clock 00:01:32)

Out[9]:

EcmResults(2RC, 3 SoC x 2 T x 2 dch / 1 chg rates, aged (SoH 92.9 % at cycle 200))

degradation_state records the aged state the fit ran at.

In [10]:

aged.degradation_state

aged.degradation_state

Out[10]:

{'lli': 0.06782521421896859,
 'lam_ne': 0.03249433822089409,
 'lam_pe': 0.0,
 'soh_pct': 92.93400366763026,
 'at_cycle': 200}

3. Fresh against aged¶

The aged grid points all sit inside the fresh grid, so the tables compare point for point wherever both fits have one.

In [11]:

temp_index = 1  # 25 degC in the grid above
rate_index = 0  # 1C discharge pulses

fig = go.Figure()
for label, ecm in (("fresh", fresh), ("aged", aged)):
    fig.add_trace(
        go.Scatter(
            x=ecm.soc_grid,
            y=1000 * ecm.r0("discharge")[:, temp_index, rate_index],
            mode="lines+markers",
            name=label,
        )
    )
fig.update_layout(
    title="R0 vs SoC at 25 degC, 1C discharge pulses",
    xaxis_title="SoC",
    yaxis_title="R0 [mOhm]",
    legend_title="Cell state",
)
fig

temp_index = 1 # 25 degC in the grid above rate_index = 0 # 1C discharge pulses fig = go.Figure() for label, ecm in (("fresh", fresh), ("aged", aged)): fig.add_trace( go.Scatter( x=ecm.soc_grid, y=1000 * ecm.r0("discharge")[:, temp_index, rate_index], mode="lines+markers", name=label, ) ) fig.update_layout( title="R0 vs SoC at 25 degC, 1C discharge pulses", xaxis_title="SoC", yaxis_title="R0 [mOhm]", legend_title="Cell state", ) fig

A small helper rebuilds the pulse voltage from the fitted tables. Feeding the same current pulse through both fits shows what the resistance growth feels like at the terminals.

In [12]:

import numpy as np


def pulse_voltage(
    ecm, soc_index, temp_index, rate_index, amps=None, t_pulse_s=60.0, t_rest_s=180.0
):
    """Modelled terminal voltage for one discharge pulse followed by a rest."""
    if amps is None:
        amps = ecm.discharge_c_rates[rate_index] * ecm.capacity
    v_ocv = float(
        np.interp(ecm.soc_grid[soc_index], ecm.ocv["soc"], ecm.ocv["mean_voltage_V"])
    )
    time_s = np.linspace(0.0, t_pulse_s + t_rest_s, 600)
    on = time_s <= t_pulse_s
    voltage = np.full_like(time_s, v_ocv)
    voltage[on] -= amps * ecm.r0("discharge")[soc_index, temp_index, rate_index]
    for k in range(1, ecm.n_rc + 1):
        r_k = ecm.r(k, "discharge")[soc_index, temp_index, rate_index]
        tau_k = r_k * ecm.c(k, "discharge")[soc_index, temp_index, rate_index]
        voltage[on] -= amps * r_k * (1.0 - np.exp(-time_s[on] / tau_k))
        u_end = amps * r_k * (1.0 - np.exp(-t_pulse_s / tau_k))
        voltage[~on] -= u_end * np.exp(-(time_s[~on] - t_pulse_s) / tau_k)
    return time_s, voltage, v_ocv

import numpy as np def pulse_voltage( ecm, soc_index, temp_index, rate_index, amps=None, t_pulse_s=60.0, t_rest_s=180.0 ): """Modelled terminal voltage for one discharge pulse followed by a rest.""" if amps is None: amps = ecm.discharge_c_rates[rate_index] * ecm.capacity v_ocv = float( np.interp(ecm.soc_grid[soc_index], ecm.ocv["soc"], ecm.ocv["mean_voltage_V"]) ) time_s = np.linspace(0.0, t_pulse_s + t_rest_s, 600) on = time_s <= t_pulse_s voltage = np.full_like(time_s, v_ocv) voltage[on] -= amps * ecm.r0("discharge")[soc_index, temp_index, rate_index] for k in range(1, ecm.n_rc + 1): r_k = ecm.r(k, "discharge")[soc_index, temp_index, rate_index] tau_k = r_k * ecm.c(k, "discharge")[soc_index, temp_index, rate_index] voltage[on] -= amps * r_k * (1.0 - np.exp(-time_s[on] / tau_k)) u_end = amps * r_k * (1.0 - np.exp(-t_pulse_s / tau_k)) voltage[~on] -= u_end * np.exp(-(time_s[~on] - t_pulse_s) / tau_k) return time_s, voltage, v_ocv

In [13]:

amps = fresh.discharge_c_rates[rate_index] * fresh.capacity

fig = go.Figure()
for label, ecm in (("fresh", fresh), ("aged", aged)):
    soc_index = int(np.argmin(np.abs(ecm.soc_grid - 0.5)))  # 0.5 SoC in each grid
    time_s, voltage, _ = pulse_voltage(
        ecm, soc_index, temp_index, rate_index, amps=amps
    )
    fig.add_trace(go.Scatter(x=time_s, y=voltage, mode="lines", name=label))
fig.update_layout(
    title="1C discharge pulse at 0.5 SoC and 25 degC, fresh against aged",
    xaxis_title="Time [s]",
    yaxis_title="Voltage [V]",
    legend_title="Cell state",
)
fig

amps = fresh.discharge_c_rates[rate_index] * fresh.capacity fig = go.Figure() for label, ecm in (("fresh", fresh), ("aged", aged)): soc_index = int(np.argmin(np.abs(ecm.soc_grid - 0.5))) # 0.5 SoC in each grid time_s, voltage, _ = pulse_voltage( ecm, soc_index, temp_index, rate_index, amps=amps ) fig.add_trace(go.Scatter(x=time_s, y=voltage, mode="lines", name=label)) fig.update_layout( title="1C discharge pulse at 0.5 SoC and 25 degC, fresh against aged", xaxis_title="Time [s]", yaxis_title="Voltage [V]", legend_title="Cell state", ) fig

The long tables take the comparison further, for example the relative R0 growth at every shared operating point.

In [14]:

fresh_table = fresh.to_dataframe().set_index(
    ["soc", "temperature_degC", "c_rate", "direction"]
)
aged_table = aged.to_dataframe().set_index(
    ["soc", "temperature_degC", "c_rate", "direction"]
)

growth = (aged_table["R0_Ohm"] / fresh_table["R0_Ohm"] - 1) * 100
growth_table = growth.dropna().rename("R0 growth [%]").reset_index()
growth_table.head(12)

fresh_table = fresh.to_dataframe().set_index( ["soc", "temperature_degC", "c_rate", "direction"] ) aged_table = aged.to_dataframe().set_index( ["soc", "temperature_degC", "c_rate", "direction"] ) growth = (aged_table["R0_Ohm"] / fresh_table["R0_Ohm"] - 1) * 100 growth_table = growth.dropna().rename("R0 growth [%]").reset_index() growth_table.head(12)

Out[14]:

	soc	temperature_degC	c_rate	direction	R0 growth [%]
0	0.1	10.0	1.0	charge	0.563724
1	0.1	10.0	1.0	discharge	0.768948
2	0.1	10.0	3.0	discharge	0.768948
3	0.1	25.0	1.0	charge	0.460181
4	0.1	25.0	1.0	discharge	0.693459
5	0.1	25.0	3.0	discharge	0.537313
6	0.5	10.0	1.0	charge	1.296899
7	0.5	10.0	1.0	discharge	1.228423
8	0.5	10.0	3.0	discharge	1.362305
9	0.5	25.0	1.0	charge	1.340987
10	0.5	25.0	1.0	discharge	1.278956
11	0.5	25.0	3.0	discharge	1.251958

The same growth as a chart, grouped by temperature.

In [15]:

at_1c = growth_table[
    (growth_table["direction"] == "discharge") & (growth_table["c_rate"] == 1.0)
]

fig = go.Figure()
for temperature, group in at_1c.groupby("temperature_degC"):
    fig.add_trace(
        go.Bar(
            x=group["soc"].astype(str),
            y=group["R0 growth [%]"],
            name=f"{temperature:g} degC",
        )
    )
fig.update_layout(
    barmode="group",
    title="R0 growth after the campaign, 1C discharge pulses",
    xaxis_title="SoC",
    yaxis_title="R0 growth [%]",
    legend_title="Temperature",
)
fig

at_1c = growth_table[ (growth_table["direction"] == "discharge") & (growth_table["c_rate"] == 1.0) ] fig = go.Figure() for temperature, group in at_1c.groupby("temperature_degC"): fig.add_trace( go.Bar( x=group["soc"].astype(str), y=group["R0 growth [%]"], name=f"{temperature:g} degC", ) ) fig.update_layout( barmode="group", title="R0 growth after the campaign, 1C discharge pulses", xaxis_title="SoC", yaxis_title="R0 growth [%]", legend_title="Temperature", ) fig

Version fe1f5cb5