Real Battery Constraints

The optimizer from Lesson 2 works, but it ignores several realities of physical batteries. If you ran that schedule on real hardware, you'd either violate electrical limits or destroy your battery prematurely. This lesson adds the constraints that make the model production-ready.

Problem: Simultaneous charge and discharge

Our Lesson 2 model allows p_charge[t] > 0 and p_discharge[t] > 0 in the same hour. Physically, a battery inverter can only flow power in one direction at a time. We need a constraint that says: in any given hour, either charge OR discharge, never both.

This requires a binary variable — turning our LP into a Mixed-Integer Program (MIP):

# Binary: 1 = discharging allowed, 0 = charging allowed
model.y = pyo.Var(model.T, within=pyo.Binary)

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

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

Minimum state of charge

Most battery management systems enforce a minimum SoC (e.g., 10%) to protect cell health and maintain reserve capacity. Similarly, many operators keep a buffer for ancillary services.

SoC_min = 0.1 * E_max  # 10% minimum
SoC_max = 0.9 * E_max  # 90% maximum — avoid stress at extremes

# Replace the simple bounds with:
model.soc = pyo.Var(model.T, bounds=(SoC_min, SoC_max))

Ramp rate limits

Batteries can ramp fast, but inverters and grid connections may limit how quickly power can change between hours. A ramp constraint limits the change in net power output:

ramp_max = 5.0  # MW per hour — max change in net output

def ramp_up_rule(m, t):
    if t == 0:
        return pyo.Constraint.Skip
    net_now = m.p_discharge[t] - m.p_charge[t]
    net_prev = m.p_discharge[t-1] - m.p_charge[t-1]
    return net_now - net_prev <= ramp_max
model.ramp_up = pyo.Constraint(model.T, rule=ramp_up_rule)

def ramp_down_rule(m, t):
    if t == 0:
        return pyo.Constraint.Skip
    net_now = m.p_discharge[t] - m.p_charge[t]
    net_prev = m.p_discharge[t-1] - m.p_charge[t-1]
    return net_prev - net_now <= ramp_max
model.ramp_down = pyo.Constraint(model.T, rule=ramp_down_rule)

Cycle-based degradation cost

Every charge/discharge cycle wears out the battery. A simple but effective model: assign a degradation cost per MWh throughput and subtract it from revenue in the objective.

# Degradation cost: €5/MWh throughput (simplified)
# Real values: derive from warranty terms or cell aging curves
deg_cost = 5.0  # €/MWh

def revenue_with_degradation(m):
    revenue = sum(prices[t] * (m.p_discharge[t] - m.p_charge[t]) for t in m.T)
    degradation = sum(deg_cost * (m.p_charge[t] + m.p_discharge[t]) for t in m.T)
    return revenue - degradation

model.objective = pyo.Objective(rule=revenue_with_degradation, sense=pyo.maximize)

End-of-horizon SoC constraint

Without a constraint, the optimizer will empty the battery by hour 24 — it has no incentive to save energy for tomorrow. Fix this by requiring SoC to return to its initial value:

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

Alternatively, add a terminal value: price the remaining energy at tomorrow's expected average price. This is more sophisticated and used in rolling-horizon optimizers.

The complete production model

"""
battery_optimizer_v2.py — With real battery constraints
"""
import pyomo.environ as pyo
import numpy as np
import pandas as pd

# --- Parameters ---
T = 24
prices = [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
ramp_max = 5.0           # MW/h
deg_cost = 5.0           # €/MWh throughput

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

# Variables
model.p_charge = pyo.Var(model.T, bounds=(0, P_max))
model.p_discharge = 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=discharge, 0=charge

# Objective: revenue minus degradation
def obj_rule(m):
    revenue = sum(prices[t] * (m.p_discharge[t] - m.p_charge[t]) for t in m.T)
    degradation = sum(deg_cost * (m.p_charge[t] + m.p_discharge[t]) for t in m.T)
    return revenue - degradation
model.objective = pyo.Objective(rule=obj_rule, sense=pyo.maximize)

# SoC balance
def soc_rule(m, t):
    soc_prev = soc_init if t == 0 else m.soc[t-1]
    return m.soc[t] == soc_prev + eta_sqrt * m.p_charge[t] - m.p_discharge[t] / eta_sqrt
model.soc_balance = pyo.Constraint(model.T, rule=soc_rule)

# No simultaneous charge/discharge
def no_simul_c(m, t):
    return m.p_charge[t] <= P_max * (1 - m.y[t])
def no_simul_d(m, t):
    return m.p_discharge[t] <= P_max * m.y[t]
model.no_simul_c = pyo.Constraint(model.T, rule=no_simul_c)
model.no_simul_d = pyo.Constraint(model.T, rule=no_simul_d)

# Ramp limits
def ramp_up(m, t):
    if t == 0: return pyo.Constraint.Skip
    return (m.p_discharge[t] - m.p_charge[t]) - (m.p_discharge[t-1] - m.p_charge[t-1]) <= ramp_max
def ramp_down(m, t):
    if t == 0: return pyo.Constraint.Skip
    return (m.p_discharge[t-1] - m.p_charge[t-1]) - (m.p_discharge[t] - m.p_charge[t]) <= ramp_max
model.ramp_up = pyo.Constraint(model.T, rule=ramp_up)
model.ramp_down = pyo.Constraint(model.T, rule=ramp_down)

# End-of-horizon: return to initial SoC
def end_soc(m):
    return m.soc[T - 1] >= soc_init
model.end_soc = pyo.Constraint(rule=end_soc)

# --- Solve ---
solver = pyo.SolverFactory('highs')
result = solver.solve(model, tee=False)

print(f"Status: {result.solver.status}")
print(f"Revenue (net of degradation): €{pyo.value(model.objective):.2f}\n")

schedule = pd.DataFrame({
    'hour': range(T),
    'price': prices,
    'charge_mw': [round(pyo.value(model.p_charge[t]), 2) for t in range(T)],
    'discharge_mw': [round(pyo.value(model.p_discharge[t]), 2) for t in range(T)],
    'soc_mwh': [round(pyo.value(model.soc[t]), 2) for t in range(T)],
    'mode': ['discharge' if pyo.value(model.y[t]) > 0.5 else 'charge/idle' for t in range(T)],
})
print(schedule.to_string(index=False))

Check your understanding

Key takeaway

Rechtliches: Impressum