simses tutorial¶

Build a BESS simulation from the simses primitives, one layer at a time:

1. Battery only — cycle with constant power, observe internal states and degradation.
2. Battery + Converter — run a small peak-shaving simulation with a Sinamics S120 loss curve.
3. Two strings — add a power-distribution rule.
4. AmbientThermalModel — couple the batteries to an ambient environment.

In [1]:

from datetime import timedelta

import numpy as np
import pandas as pd
from tqdm.notebook import tqdm

from simses.battery.battery import Battery
from simses.converter.converter import Converter
from simses.degradation.state import DegradationState
from simses.model.cell.sony_lfp import SonyLFP
from simses.model.converter.sinamics import SinamicsS120Fit
from simses.thermal.ambient import AmbientThermalModel

from datetime import timedelta import numpy as np import pandas as pd from tqdm.notebook import tqdm from simses.battery.battery import Battery from simses.converter.converter import Converter from simses.degradation.state import DegradationState from simses.model.cell.sony_lfp import SonyLFP from simses.model.converter.sinamics import SinamicsS120Fit from simses.thermal.ambient import AmbientThermalModel

1. Battery¶

A Battery wraps a CellType in a series–parallel circuit and steps forward via step(power_setpoint, dt).

We use the SonyLFP cell with its default calendar + cyclic degradation model so soh_Q / soh_R drift is visible.

In [2]:

cell = SonyLFP()

battery = Battery(
    cell=cell,
    circuit=(104, 10),  # 104 cells in series, 10 parallel strings
    initial_states={"start_soc": 0.5, "start_T": 25.0},
    degradation=SonyLFP.default_degradation_model(
        initial_soc=0.5,
        initial_state=DegradationState(qloss_cal=1e-4),
    ),
)

pd.Series({
    "nominal_capacity [Ah]": battery.nominal_capacity,
    "nominal_voltage [V]": battery.nominal_voltage,
    "nominal_energy [kWh]": battery.nominal_energy_capacity / 1e3,
    "max_charge_current [A]": battery.max_charge_current,
    "max_discharge_current [A]": battery.max_discharge_current,
    "thermal_capacity [kJ/K]": battery.thermal_capacity / 1e3,
})

cell = SonyLFP() battery = Battery( cell=cell, circuit=(104, 10), # 104 cells in series, 10 parallel strings initial_states={"start_soc": 0.5, "start_T": 25.0}, degradation=SonyLFP.default_degradation_model( initial_soc=0.5, initial_state=DegradationState(qloss_cal=1e-4), ), ) pd.Series({ "nominal_capacity [Ah]": battery.nominal_capacity, "nominal_voltage [V]": battery.nominal_voltage, "nominal_energy [kWh]": battery.nominal_energy_capacity / 1e3, "max_charge_current [A]": battery.max_charge_current, "max_discharge_current [A]": battery.max_discharge_current, "thermal_capacity [kJ/K]": battery.thermal_capacity / 1e3, })

Out[2]:

nominal_capacity [Ah]         30.0000
nominal_voltage [V]          332.8000
nominal_energy [kWh]           9.9840
max_charge_current [A]        30.0000
max_discharge_current [A]    198.0000
thermal_capacity [kJ/K]       72.8728
dtype: float64

Cycle the battery with a ±5 kW square wave (30-min half-period) for 48 h. battery.state is a plain dataclass — we log whichever fields interest us.

In [3]:

dt = 60.0
duration = timedelta(hours=48)
n_steps = int(duration.total_seconds() / dt)

half_period = 30 * 60  # 30 minutes
power_amp = 5e3
power_profile = np.where((np.arange(n_steps) * dt) % (2 * half_period) < half_period, power_amp, -power_amp)

log = {k: np.empty(n_steps) for k in ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"]}
for i, p in enumerate(tqdm(power_profile)):
    battery.step(float(p), dt)
    for k in log:
        log[k][i] = getattr(battery.state, k)

index = pd.date_range("2025-01-01", periods=n_steps, freq=f"{int(dt)}s")
df_bat = pd.DataFrame(log, index=index)
df_bat.head()

dt = 60.0 duration = timedelta(hours=48) n_steps = int(duration.total_seconds() / dt) half_period = 30 * 60 # 30 minutes power_amp = 5e3 power_profile = np.where((np.arange(n_steps) * dt) % (2 * half_period) < half_period, power_amp, -power_amp) log = {k: np.empty(n_steps) for k in ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"]} for i, p in enumerate(tqdm(power_profile)): battery.step(float(p), dt) for k in log: log[k][i] = getattr(battery.state, k) index = pd.date_range("2025-01-01", periods=n_steps, freq=f"{int(dt)}s") df_bat = pd.DataFrame(log, index=index) df_bat.head()

Out[3]:

	soc	v	i	T	loss	heat	soh_Q	soh_R	power
2025-01-01 00:00:00	0.507898	351.710099	14.216254	25.0	121.757857	170.366187	0.999984	1.0	5000.0
2025-01-01 00:01:00	0.515794	351.776486	14.213571	25.0	122.134143	171.373127	0.999970	1.0	5000.0
2025-01-01 00:02:00	0.523690	351.832855	14.211294	25.0	122.545115	172.415825	0.999957	1.0	5000.0
2025-01-01 00:03:00	0.531584	351.894555	14.208802	25.0	123.027025	173.528400	0.999946	1.0	5000.0
2025-01-01 00:04:00	0.539477	351.967217	14.205868	25.0	123.553128	174.683251	0.999935	1.0	5000.0

In [4]:

df_bat["soc"].plot(title="SOC", figsize=(10, 3));

df_bat["soc"].plot(title="SOC", figsize=(10, 3));

In [5]:

df_bat[["v"]].plot(title="Terminal voltage [V]", figsize=(10, 3));

df_bat[["v"]].plot(title="Terminal voltage [V]", figsize=(10, 3));

In [6]:

ax = df_bat[["soh_Q", "soh_R"]].plot(secondary_y="soh_R", title="State of health", figsize=(10, 3))
ax.set_ylabel("soh_Q [p.u.]")
ax.right_ax.set_ylabel("soh_R [p.u.]");

ax = df_bat[["soh_Q", "soh_R"]].plot(secondary_y="soh_R", title="State of health", figsize=(10, 3)) ax.set_ylabel("soh_Q [p.u.]") ax.right_ax.set_ylabel("soh_R [p.u.]");

In [7]:

df_bat[["loss", "heat"]].plot(title="Battery losses and heat [W]", figsize=(10, 3));

df_bat[["loss", "heat"]].plot(title="Battery losses and heat [W]", figsize=(10, 3));

2. Battery + Converter¶

Converter takes a downstream storage (here: a Battery), a rated max_power, and a loss_model. SinamicsS120Fit is a parametric fit of the Siemens S120 efficiency curve — nonlinear around the low-power region.

In [8]:

battery_1 = Battery(
    cell=cell,
    circuit=(104, 10),
    initial_states={"start_soc": 0.5, "start_T": 25.0},
    degradation=True,
)

converter_1 = Converter(
    loss_model=SinamicsS120Fit(),
    max_power=20e3,  # 20 kW
    storage=battery_1,
)

battery_1 = Battery( cell=cell, circuit=(104, 10), initial_states={"start_soc": 0.5, "start_T": 25.0}, degradation=True, ) converter_1 = Converter( loss_model=SinamicsS120Fit(), max_power=20e3, # 20 kW storage=battery_1, )

Sweep AC power across [-max_power, +max_power] and look at the efficiency curve the converter sees. Note the efficiency dip at very low power.

In [9]:

p_ac = np.linspace(-converter_1.max_power, converter_1.max_power, 401)
p_dc = np.array([converter_1.ac_to_dc(p) for p in p_ac])

# efficiency = |useful| / |input|; charging -> dc useful, discharging -> ac useful
eta = np.where(p_ac > 0, np.abs(p_dc) / np.abs(p_ac), np.abs(p_ac) / np.abs(p_dc))

df_eff = pd.DataFrame({"efficiency": eta, "loss [W]": p_ac - p_dc}, index=pd.Index(p_ac / 1e3, name="AC power [kW]"))
df_eff.plot(subplots=True, figsize=(10, 5), title="SinamicsS120Fit: efficiency and loss vs AC power");

p_ac = np.linspace(-converter_1.max_power, converter_1.max_power, 401) p_dc = np.array([converter_1.ac_to_dc(p) for p in p_ac]) # efficiency = |useful| / |input|; charging -> dc useful, discharging -> ac useful eta = np.where(p_ac > 0, np.abs(p_dc) / np.abs(p_ac), np.abs(p_ac) / np.abs(p_dc)) df_eff = pd.DataFrame({"efficiency": eta, "loss [W]": p_ac - p_dc}, index=pd.Index(p_ac / 1e3, name="AC power [kW]")) df_eff.plot(subplots=True, figsize=(10, 5), title="SinamicsS120Fit: efficiency and loss vs AC power");

/tmp/ipykernel_2435/686983845.py:5: RuntimeWarning: invalid value encountered in divide
  eta = np.where(p_ac > 0, np.abs(p_dc) / np.abs(p_ac), np.abs(p_ac) / np.abs(p_dc))

Build a synthetic 48 h load profile (daily pattern + peaks + noise), scaled to 20 kW peak. Then run a trivial peak-shaving EMS at a 15 kW threshold.

In [10]:

class SimplePeakShaving:
    """Discharge when load exceeds threshold, charge otherwise (until nearly full)."""

    def __init__(self, power_threshold: float, soc_threshold: float = 0.999):
        self.power_threshold = power_threshold
        self.soc_threshold = soc_threshold

    def next(self, power: float, soc: float) -> float:
        net_power = self.power_threshold - power
        if net_power < 0.0:
            return net_power  # discharge
        if soc < self.soc_threshold:
            return net_power  # charge
        return 0.0

class SimplePeakShaving: """Discharge when load exceeds threshold, charge otherwise (until nearly full).""" def __init__(self, power_threshold: float, soc_threshold: float = 0.999): self.power_threshold = power_threshold self.soc_threshold = soc_threshold def next(self, power: float, soc: float) -> float: net_power = self.power_threshold - power if net_power < 0.0: return net_power # discharge if soc < self.soc_threshold: return net_power # charge return 0.0

In [11]:

def make_load_profile(hours=48, dt=60.0, peak_power=20e3, seed=0):
    rng = np.random.default_rng(seed)
    n = int(hours * 3600 / dt)
    t = np.arange(n) * dt
    daily = 0.45 + 0.35 * np.sin(2 * np.pi * (t / 86400 - 0.25))  # low at night, high midday
    peaks = np.zeros(n)
    for _ in range(hours // 4):  # a handful of demand spikes
        start = rng.integers(0, n - int(1800 / dt))
        width = int(rng.integers(300, 1800) / dt)
        peaks[start : start + width] += rng.uniform(0.3, 0.6)
    noise = 0.03 * rng.standard_normal(n)
    profile = np.clip(daily + peaks + noise, 0.05, None)
    profile *= peak_power / profile.max()
    return pd.DataFrame({"load": profile}, index=pd.date_range("2025-01-01", periods=n, freq=f"{int(dt)}s"))


load_profile = make_load_profile(hours=48, dt=dt, peak_power=20e3)
load_profile["load"].plot(title="Synthetic load [W]", figsize=(10, 3));

def make_load_profile(hours=48, dt=60.0, peak_power=20e3, seed=0): rng = np.random.default_rng(seed) n = int(hours * 3600 / dt) t = np.arange(n) * dt daily = 0.45 + 0.35 * np.sin(2 * np.pi * (t / 86400 - 0.25)) # low at night, high midday peaks = np.zeros(n) for _ in range(hours // 4): # a handful of demand spikes start = rng.integers(0, n - int(1800 / dt)) width = int(rng.integers(300, 1800) / dt) peaks[start : start + width] += rng.uniform(0.3, 0.6) noise = 0.03 * rng.standard_normal(n) profile = np.clip(daily + peaks + noise, 0.05, None) profile *= peak_power / profile.max() return pd.DataFrame({"load": profile}, index=pd.date_range("2025-01-01", periods=n, freq=f"{int(dt)}s")) load_profile = make_load_profile(hours=48, dt=dt, peak_power=20e3) load_profile["load"].plot(title="Synthetic load [W]", figsize=(10, 3));

In [12]:

ems = SimplePeakShaving(power_threshold=15e3)

log = {k: np.empty(len(load_profile)) for k in ["load", "setpoint", "ac_power", "soc", "conv_loss", "bat_loss"]}
for i, load in enumerate(tqdm(load_profile["load"].to_numpy())):
    setpoint = ems.next(load, battery_1.state.soc)
    converter_1.step(setpoint, dt)
    log["load"][i] = load
    log["setpoint"][i] = setpoint
    log["ac_power"][i] = converter_1.state.power
    log["soc"][i] = battery_1.state.soc
    log["conv_loss"][i] = converter_1.state.loss
    log["bat_loss"][i] = battery_1.state.loss

df_ps = pd.DataFrame(log, index=load_profile.index)
df_ps["grid"] = df_ps["load"] + df_ps["ac_power"]

ems = SimplePeakShaving(power_threshold=15e3) log = {k: np.empty(len(load_profile)) for k in ["load", "setpoint", "ac_power", "soc", "conv_loss", "bat_loss"]} for i, load in enumerate(tqdm(load_profile["load"].to_numpy())): setpoint = ems.next(load, battery_1.state.soc) converter_1.step(setpoint, dt) log["load"][i] = load log["setpoint"][i] = setpoint log["ac_power"][i] = converter_1.state.power log["soc"][i] = battery_1.state.soc log["conv_loss"][i] = converter_1.state.loss log["bat_loss"][i] = battery_1.state.loss df_ps = pd.DataFrame(log, index=load_profile.index) df_ps["grid"] = df_ps["load"] + df_ps["ac_power"]

In [13]:

ax = df_ps[["load", "grid"]].plot(title="Peak shaving at 15 kW", figsize=(10, 3))
ax.axhline(ems.power_threshold, color="red", linestyle="--", linewidth=0.8);

ax = df_ps[["load", "grid"]].plot(title="Peak shaving at 15 kW", figsize=(10, 3)) ax.axhline(ems.power_threshold, color="red", linestyle="--", linewidth=0.8);

In [14]:

df_ps["soc"].plot(title="SOC", figsize=(10, 3));

df_ps["soc"].plot(title="SOC", figsize=(10, 3));

In [15]:

df_ps[["conv_loss", "bat_loss"]].plot(title="Converter vs battery losses [W]", figsize=(10, 3));

df_ps[["conv_loss", "bat_loss"]].plot(title="Converter vs battery losses [W]", figsize=(10, 3));

3. Two strings with power distribution¶

Each string is an independent (Battery, Converter) pair. To make the split non-trivial, the second string has half the capacity and half the converter rating of the first. The strings also start at different SOCs (0.3 and 0.7) so the SOC-weighted distribution rule can visibly drain the full string first and fill the empty one first.

In [16]:

def make_string(start_soc, circuit=(104, 10), max_power=20e3, start_T=25.0, effective_cooling_area=1.0):
    bat = Battery(
        cell=cell,
        circuit=circuit,
        initial_states={"start_soc": start_soc, "start_T": start_T},
        degradation=True,
        effective_cooling_area=effective_cooling_area,
    )
    conv = Converter(loss_model=SinamicsS120Fit(), max_power=max_power, storage=bat)
    return bat, conv


# s0: full-size string; s1: half capacity, half power
strings = {
    "s0": make_string(start_soc=0.3, circuit=(104, 10), max_power=20e3),
    "s1": make_string(start_soc=0.7, circuit=(104, 5), max_power=10e3),
}

def make_string(start_soc, circuit=(104, 10), max_power=20e3, start_T=25.0, effective_cooling_area=1.0): bat = Battery( cell=cell, circuit=circuit, initial_states={"start_soc": start_soc, "start_T": start_T}, degradation=True, effective_cooling_area=effective_cooling_area, ) conv = Converter(loss_model=SinamicsS120Fit(), max_power=max_power, storage=bat) return bat, conv # s0: full-size string; s1: half capacity, half power strings = { "s0": make_string(start_soc=0.3, circuit=(104, 10), max_power=20e3), "s1": make_string(start_soc=0.7, circuit=(104, 5), max_power=10e3), }

In [17]:

def power_distribution(strings: dict, target_power: float) -> dict:
    """SOC-weighted split: discharge drains the full string first, charge fills the empty one first."""
    eps = 1e-6
    socs = np.array([bat.state.soc for bat, _ in strings.values()])
    if target_power < 0:  # discharging -> prioritise high SOC
        weights = (socs + eps) / np.sum(socs + eps)
    else:  # charging -> prioritise low SOC
        weights = (1 - socs + eps) / np.sum(1 - socs + eps)
    return {string: target_power * weights[i] for i, string in enumerate(strings)}

def power_distribution(strings: dict, target_power: float) -> dict: """SOC-weighted split: discharge drains the full string first, charge fills the empty one first.""" eps = 1e-6 socs = np.array([bat.state.soc for bat, _ in strings.values()]) if target_power < 0: # discharging -> prioritise high SOC weights = (socs + eps) / np.sum(socs + eps) else: # charging -> prioritise low SOC weights = (1 - socs + eps) / np.sum(1 - socs + eps) return {string: target_power * weights[i] for i, string in enumerate(strings)}

In [18]:

ems = SimplePeakShaving(power_threshold=30e3)
load_profile_2 = make_load_profile(hours=48, dt=dt, peak_power=40e3, seed=1)

cols = ["load", "setpoint", "grid_power"]
for key in strings:
    cols += [f"{key}_soc", f"{key}_power", f"{key}_T"]
log = {k: np.empty(len(load_profile_2)) for k in cols}

for i, load in enumerate(tqdm(load_profile_2["load"].to_numpy())):
    avg_soc = float(np.mean([bat.state.soc for bat, _ in strings.values()]))
    setpoint = ems.next(load, avg_soc)

    powers = power_distribution(strings, setpoint)
    total_ac = 0.0
    for key, (bat, conv) in strings.items():
        conv.step(powers[key], dt)
        log[f"{key}_soc"][i] = bat.state.soc
        log[f"{key}_power"][i] = conv.state.power
        log[f"{key}_T"][i] = bat.state.T
        total_ac += conv.state.power

    log["load"][i] = load
    log["setpoint"][i] = setpoint
    log["grid_power"][i] = load + total_ac

df2 = pd.DataFrame(log, index=load_profile_2.index)

ems = SimplePeakShaving(power_threshold=30e3) load_profile_2 = make_load_profile(hours=48, dt=dt, peak_power=40e3, seed=1) cols = ["load", "setpoint", "grid_power"] for key in strings: cols += [f"{key}_soc", f"{key}_power", f"{key}_T"] log = {k: np.empty(len(load_profile_2)) for k in cols} for i, load in enumerate(tqdm(load_profile_2["load"].to_numpy())): avg_soc = float(np.mean([bat.state.soc for bat, _ in strings.values()])) setpoint = ems.next(load, avg_soc) powers = power_distribution(strings, setpoint) total_ac = 0.0 for key, (bat, conv) in strings.items(): conv.step(powers[key], dt) log[f"{key}_soc"][i] = bat.state.soc log[f"{key}_power"][i] = conv.state.power log[f"{key}_T"][i] = bat.state.T total_ac += conv.state.power log["load"][i] = load log["setpoint"][i] = setpoint log["grid_power"][i] = load + total_ac df2 = pd.DataFrame(log, index=load_profile_2.index)

In [19]:

df2[["s0_soc", "s1_soc"]].plot(title="SOC of both strings", figsize=(10, 3));

df2[["s0_soc", "s1_soc"]].plot(title="SOC of both strings", figsize=(10, 3));

In [20]:

df2[["s0_power", "s1_power"]].plot(title="Per-string AC power [W]", figsize=(10, 3));

df2[["s0_power", "s1_power"]].plot(title="Per-string AC power [W]", figsize=(10, 3));

In [21]:

ax = df2[["load", "grid_power"]].plot(title="Peak shaving at 30 kW (two strings)", figsize=(10, 3))
ax.axhline(30e3, color="red", linestyle="--", linewidth=0.8);

ax = df2[["load", "grid_power"]].plot(title="Peak shaving at 30 kW (two strings)", figsize=(10, 3)) ax.axhline(30e3, color="red", linestyle="--", linewidth=0.8);

4. Thermal model¶

AmbientThermalModel is the simplest thermal environment: each registered component exchanges heat with a constant-temperature ambient. Registration is pure duck typing — any object exposing state.T, state.heat, thermal_capacity, thermal_resistance works.

To make the thermal dynamics visible we restrict each battery's effective_cooling_area to 10% of the cell surface — roughly what you get with a single-sided cooling plate.

For richer setups (container + HVAC, room with solar gains) see ContainerThermalModel.

In [22]:

strings_th = {
    "s0": make_string(start_soc=0.3, circuit=(104, 10), max_power=20e3, effective_cooling_area=0.1),
    "s1": make_string(start_soc=0.7, circuit=(104, 5), max_power=10e3, effective_cooling_area=0.1),
}

thermal = AmbientThermalModel(T_ambient=25.0)
for bat, _ in strings_th.values():
    thermal.add_component(bat)

strings_th = { "s0": make_string(start_soc=0.3, circuit=(104, 10), max_power=20e3, effective_cooling_area=0.1), "s1": make_string(start_soc=0.7, circuit=(104, 5), max_power=10e3, effective_cooling_area=0.1), } thermal = AmbientThermalModel(T_ambient=25.0) for bat, _ in strings_th.values(): thermal.add_component(bat)

In [23]:

ems = SimplePeakShaving(power_threshold=30e3)

cols = ["load"]
for key in strings_th:
    cols += [f"{key}_soc", f"{key}_T", f"{key}_heat", f"{key}_soh_Q"]
log = {k: np.empty(len(load_profile_2)) for k in cols}

for i, load in enumerate(tqdm(load_profile_2["load"].to_numpy())):
    avg_soc = float(np.mean([bat.state.soc for bat, _ in strings_th.values()]))
    setpoint = ems.next(load, avg_soc)

    powers = power_distribution(strings_th, setpoint)
    for key, (_bat, conv) in strings_th.items():
        conv.step(powers[key], dt)
    thermal.step(dt)  # thermal step *after* the electrical step

    for key, (bat, _) in strings_th.items():
        log[f"{key}_soc"][i] = bat.state.soc
        log[f"{key}_T"][i] = bat.state.T
        log[f"{key}_heat"][i] = bat.state.heat
        log[f"{key}_soh_Q"][i] = bat.state.soh_Q
    log["load"][i] = load

df3 = pd.DataFrame(log, index=load_profile_2.index)

ems = SimplePeakShaving(power_threshold=30e3) cols = ["load"] for key in strings_th: cols += [f"{key}_soc", f"{key}_T", f"{key}_heat", f"{key}_soh_Q"] log = {k: np.empty(len(load_profile_2)) for k in cols} for i, load in enumerate(tqdm(load_profile_2["load"].to_numpy())): avg_soc = float(np.mean([bat.state.soc for bat, _ in strings_th.values()])) setpoint = ems.next(load, avg_soc) powers = power_distribution(strings_th, setpoint) for key, (_bat, conv) in strings_th.items(): conv.step(powers[key], dt) thermal.step(dt) # thermal step *after* the electrical step for key, (bat, _) in strings_th.items(): log[f"{key}_soc"][i] = bat.state.soc log[f"{key}_T"][i] = bat.state.T log[f"{key}_heat"][i] = bat.state.heat log[f"{key}_soh_Q"][i] = bat.state.soh_Q log["load"][i] = load df3 = pd.DataFrame(log, index=load_profile_2.index)

In [24]:

df3[["s0_T", "s1_T"]].plot(title="Battery temperature [°C]", figsize=(10, 3));

df3[["s0_T", "s1_T"]].plot(title="Battery temperature [°C]", figsize=(10, 3));

In [25]:

df3[["s0_heat", "s1_heat"]].plot(title="Heat generation [W]", figsize=(10, 3));

df3[["s0_heat", "s1_heat"]].plot(title="Heat generation [W]", figsize=(10, 3));

In [26]:

df3[["s0_soh_Q", "s1_soh_Q"]].plot(title="Capacity SOH", figsize=(10, 3));

df3[["s0_soh_Q", "s1_soh_Q"]].plot(title="Capacity SOH", figsize=(10, 3));