Reale Batterie-Nebenbedingungen

Der Optimierer aus Lektion 2 funktioniert, ignoriert aber mehrere Realitäten physischer Batterien. Würdest du diesen Fahrplan auf echter Hardware fahren, würdest du entweder elektrische Grenzen verletzen oder die Batterie vorzeitig zerstören. Diese Lektion fügt die Nebenbedingungen hinzu, die das Modell produktionsreif machen.

Problem: Gleichzeitiges Laden und Entladen

Unser Lektion-2-Modell erlaubt p_laden[t] > 0 und p_entladen[t] > 0 in derselben Stunde. Physisch kann ein Batterie-Wechselrichter nur in eine Richtung Leistung fließen lassen. Wir brauchen eine Nebenbedingung: In jeder Stunde entweder Laden ODER Entladen, niemals beides.

Das erfordert eine Binärvariable — unser LP wird zum Mixed-Integer Program (MIP):

# Binär: 1 = Entladen erlaubt, 0 = Laden erlaubt
model.y = pyo.Var(model.T, within=pyo.Binary)

def kein_gleichz_laden(m, t):
    return m.p_laden[t] <= P_max * (1 - m.y[t])
model.kein_gleichz_laden = pyo.Constraint(model.T, rule=kein_gleichz_laden)

def kein_gleichz_entladen(m, t):
    return m.p_entladen[t] <= P_max * m.y[t]
model.kein_gleichz_entladen = pyo.Constraint(model.T, rule=kein_gleichz_entladen)

Minimaler Ladezustand

Die meisten Batteriemanagementsysteme erzwingen ein minimales SoC (z.B. 10%) zum Schutz der Zellgesundheit und zur Aufrechterhaltung einer Reservekapazität. Ebenso halten viele Betreiber einen Puffer für Regelleistung.

SoC_min = 0.1 * E_max  # 10% Minimum
SoC_max = 0.9 * E_max  # 90% Maximum — Stress an den Extremen vermeiden

# Ersetze die einfachen Grenzen mit:
model.soc = pyo.Var(model.T, bounds=(SoC_min, SoC_max))

Rampenratenbegrenzung

Batterien können schnell rampen, aber Wechselrichter und Netzanschlüsse können begrenzen, wie schnell sich die Leistung zwischen Stunden ändern darf. Eine Rampennebenbedingung begrenzt die Änderung der Nettoleistung:

rampe_max = 5.0  # MW pro Stunde — max. Änderung der Nettoleistung

def rampe_hoch_regel(m, t):
    if t == 0:
        return pyo.Constraint.Skip
    netto_jetzt = m.p_entladen[t] - m.p_laden[t]
    netto_vorher = m.p_entladen[t-1] - m.p_laden[t-1]
    return netto_jetzt - netto_vorher <= rampe_max
model.rampe_hoch = pyo.Constraint(model.T, rule=rampe_hoch_regel)

def rampe_runter_regel(m, t):
    if t == 0:
        return pyo.Constraint.Skip
    netto_jetzt = m.p_entladen[t] - m.p_laden[t]
    netto_vorher = m.p_entladen[t-1] - m.p_laden[t-1]
    return netto_vorher - netto_jetzt <= rampe_max
model.rampe_runter = pyo.Constraint(model.T, rule=rampe_runter_regel)

Zyklusbasierte Degradationskosten

Jeder Lade-/Entladezyklus verschleißt die Batterie. Ein einfaches aber effektives Modell: weise einen Degradationskostensatz pro MWh Durchsatz zu und ziehe ihn in der Zielfunktion vom Erlös ab.

# Degradationskosten: €5/MWh Durchsatz (vereinfacht)
# Reale Werte: aus Garantiebedingungen oder Zellalterungskurven ableiten
deg_kosten = 5.0  # €/MWh

def erloes_mit_degradation(m):
    erloes = sum(preise[t] * (m.p_entladen[t] - m.p_laden[t]) for t in m.T)
    degradation = sum(deg_kosten * (m.p_laden[t] + m.p_entladen[t]) for t in m.T)
    return erloes - degradation

model.zielfunktion = pyo.Objective(rule=erloes_mit_degradation, sense=pyo.maximize)

End-of-Horizon SoC-Nebenbedingung

Ohne Nebenbedingung wird der Optimierer die Batterie bis Stunde 24 leeren — er hat keinen Anreiz, Energie für morgen aufzubewahren. Behebe das, indem du verlangst, dass SoC auf den Anfangswert zurückkehrt:

def end_soc_regel(m):
    return m.soc[T - 1] >= soc_init
model.end_soc = pyo.Constraint(rule=end_soc_regel)

Alternativ füge einen Endwert hinzu: bewerte die verbleibende Energie zum erwarteten Durchschnittspreis von morgen. Das ist ausgefeilter und wird in rollierenden Optimierern verwendet.

Das vollständige Produktionsmodell

"""
batterie_optimierer_v2.py — Mit realen Batterie-Nebenbedingungen
"""
import pyomo.environ as pyo
import numpy as np
import pandas as pd

# --- Parameter ---
T = 24
preise = [28, 25, 22, 20, 19, 21, 35, 55,
          62, 58, 52, 48, 45, 44, 47, 53,
          68, 85, 92, 78, 60, 45, 35, 30]

P_max = 10.0
E_max = 20.0
eta = 0.90
eta_sqrt = np.sqrt(eta)
soc_init = 10.0
SoC_min = 0.1 * E_max   # 2 MWh
SoC_max = 0.9 * E_max   # 18 MWh
rampe_max = 5.0          # MW/h
deg_kosten = 5.0         # €/MWh Durchsatz

# --- Modell ---
model = pyo.ConcreteModel()
model.T = pyo.RangeSet(0, T - 1)

# Variablen
model.p_laden = pyo.Var(model.T, bounds=(0, P_max))
model.p_entladen = pyo.Var(model.T, bounds=(0, P_max))
model.soc = pyo.Var(model.T, bounds=(SoC_min, SoC_max))
model.y = pyo.Var(model.T, within=pyo.Binary)  # 1=Entladen, 0=Laden

# Zielfunktion: Erlös minus Degradation
def ziel_regel(m):
    erloes = sum(preise[t] * (m.p_entladen[t] - m.p_laden[t]) for t in m.T)
    degradation = sum(deg_kosten * (m.p_laden[t] + m.p_entladen[t]) for t in m.T)
    return erloes - degradation
model.zielfunktion = pyo.Objective(rule=ziel_regel, sense=pyo.maximize)

# SoC-Bilanz
def soc_regel(m, t):
    soc_vorher = soc_init if t == 0 else m.soc[t-1]
    return m.soc[t] == soc_vorher + eta_sqrt * m.p_laden[t] - m.p_entladen[t] / eta_sqrt
model.soc_bilanz = pyo.Constraint(model.T, rule=soc_regel)

# Kein gleichzeitiges Laden/Entladen
def kgl_l(m, t):
    return m.p_laden[t] <= P_max * (1 - m.y[t])
def kgl_e(m, t):
    return m.p_entladen[t] <= P_max * m.y[t]
model.kgl_l = pyo.Constraint(model.T, rule=kgl_l)
model.kgl_e = pyo.Constraint(model.T, rule=kgl_e)

# Rampengrenzen
def rampe_h(m, t):
    if t == 0: return pyo.Constraint.Skip
    return (m.p_entladen[t] - m.p_laden[t]) - (m.p_entladen[t-1] - m.p_laden[t-1]) <= rampe_max
def rampe_r(m, t):
    if t == 0: return pyo.Constraint.Skip
    return (m.p_entladen[t-1] - m.p_laden[t-1]) - (m.p_entladen[t] - m.p_laden[t]) <= rampe_max
model.rampe_h = pyo.Constraint(model.T, rule=rampe_h)
model.rampe_r = pyo.Constraint(model.T, rule=rampe_r)

# End-of-Horizon: zurück zum Anfangs-SoC
def end_soc(m):
    return m.soc[T - 1] >= soc_init
model.end_soc = pyo.Constraint(rule=end_soc)

# --- Lösen ---
solver = pyo.SolverFactory('highs')
ergebnis = solver.solve(model, tee=False)

print(f"Status: {ergebnis.solver.status}")
print(f"Erlös (netto nach Degradation): €{pyo.value(model.zielfunktion):.2f}\n")

fahrplan = pd.DataFrame({
    'stunde': range(T),
    'preis': preise,
    'laden_mw': [round(pyo.value(model.p_laden[t]), 2) for t in range(T)],
    'entladen_mw': [round(pyo.value(model.p_entladen[t]), 2) for t in range(T)],
    'soc_mwh': [round(pyo.value(model.soc[t]), 2) for t in range(T)],
    'modus': ['Entladen' if pyo.value(model.y[t]) > 0.5 else 'Laden/Leerlauf' for t in range(T)],
})
print(fahrplan.to_string(index=False))

Überprüfe dein Verständnis

Rechtliches: Impressum