pluginbattery/scripts/capacity_sweep.py
Michiel Berger 60e0706736 Initial import: home-battery ROI simulator + cracked thuisbatterijgids calc
- Hourly data exporter (InfluxDB → CSV) for prices, P1, irradiance.
- LP-based 24h-foresight oracle dispatch with SoC-consistent state engine.
- Reverse-engineered thuisbatterijgids.nl formula (matches their quotes
  to within €0.50 across three battery configs).
- Catalog scraper for the 52 batteries on thuisbatterijgids.net via their
  /wp-json REST endpoint.
- Web app (Flask) that ranks every catalog battery by honest payback and
  contrasts with the store's quote, deployable via the included Procfile.
2026-04-30 13:46:27 +02:00

126 lines
4.7 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#!/usr/bin/env python3
"""Sweep battery capacity at fixed power and report savings.
Same scenario as the dynamic-rate calibration test:
- 4000 kWh/yr demand
- Dynamic rate, avg retail €0.25
- No PV, no saldering
- 0.8 kW charge/discharge power (plug-in)
Reports both the honest LP year-1 savings and the cracked store-calculator
quote, plus the marginal savings per added kWh of capacity.
"""
from __future__ import annotations
import argparse
import csv
import math
from pathlib import Path
import numpy as np
from pluginbattery.sim import (
Battery,
apply_nl_tariff,
load_hourly,
oracle_daily_schedule,
simulate,
)
from pluginbattery.store_calc import (
Scenario,
StoreParams,
quote as store_quote,
)
def main() -> None:
p = argparse.ArgumentParser(description=__doc__)
p.add_argument("--power", type=float, default=0.8, help="Battery in/out power kW")
p.add_argument("--retail", type=float, default=0.25, help="Avg retail EUR/kWh")
p.add_argument("--demand", type=float, default=4000.0, help="Annual demand kWh/yr")
p.add_argument("--eta", type=float, default=0.88, help="LP round-trip efficiency")
p.add_argument("--cap-min", type=float, default=0.5, help="Min capacity kWh")
p.add_argument("--cap-max", type=float, default=15.0, help="Max capacity kWh")
p.add_argument("--cap-step", type=float, default=0.5, help="Capacity step kWh")
p.add_argument("--out", default="data/processed/capacity_sweep.csv")
args = p.parse_args()
base = load_hourly("data/raw"); base = apply_nl_tariff(base)
# Scale to the test scenario.
annual_demand = base["demand_kwh"].sum() * (8766.0 / len(base))
df = base.copy()
df["demand_kwh"] = base["demand_kwh"] * (args.demand / annual_demand)
scale = args.retail / base["eur_per_kwh"].mean()
df["eur_per_kwh"] = base["eur_per_kwh"] * scale
df["epex_eur_per_kwh"] = base["epex_eur_per_kwh"] * scale
df["export_eur_per_kwh"] = 0.0 # no saldering
avg_epex = float(df["epex_eur_per_kwh"].mean())
store_params = StoreParams(avg_epex_eur_per_kwh=avg_epex)
capacities = np.arange(args.cap_min, args.cap_max + args.cap_step / 2, args.cap_step)
rows = []
print(f"Sweeping {len(capacities)} capacities at {args.power} kW power...")
print()
print(f"{'cap (kWh)':>9s} {'LP €/yr':>9s} {'store €/yr':>11s} "
f"{'LP Δ/kWh':>10s} {'store Δ/kWh':>12s} {'overstate':>10s}")
print("-" * 75)
prev_lp, prev_store = 0.0, 0.0
for cap in capacities:
bat = Battery(
capacity_kwh=float(cap), max_charge_kw=args.power,
max_discharge_kw=args.power, round_trip_eff=args.eta,
allows_export=False,
)
out = simulate(df, bat, oracle_daily_schedule(df, bat))
lp = float(out["savings"].sum())
s = Scenario(
capacity_kwh=float(cap), max_charge_kw=args.power,
battery_cost_eur=0.0, # not used here
avg_retail_eur_per_kwh=args.retail,
has_pv=False, has_saldering=False, dynamic_rate=True,
)
sc = store_quote(s, store_params)["year1_total"]
marg_lp = (lp - prev_lp) / args.cap_step
marg_store = (sc - prev_store) / args.cap_step
ratio = sc / lp if lp > 0 else float("inf")
prev_lp, prev_store = lp, sc
rows.append({
"capacity_kwh": round(float(cap), 3),
"power_kw": args.power,
"lp_year1_eur": round(lp, 2),
"store_year1_eur": round(sc, 2),
"marginal_lp_eur_per_kwh": round(marg_lp, 2),
"marginal_store_eur_per_kwh": round(marg_store, 2),
"overstatement_ratio": round(ratio, 2),
})
print(f"{cap:>9.2f}{lp:>7.2f}{sc:>9.2f} "
f"{marg_lp:>8.2f}{marg_store:>10.2f} {ratio:>8.2f}×")
out_path = Path(args.out)
out_path.parent.mkdir(parents=True, exist_ok=True)
with out_path.open("w", newline="") as f:
writer = csv.DictWriter(f, fieldnames=rows[0].keys())
writer.writeheader()
writer.writerows(rows)
print(f"\n{out_path}")
# Identify the diminishing-returns elbow on the LP curve.
lp_marg = [r["marginal_lp_eur_per_kwh"] for r in rows]
threshold = lp_marg[0] * 0.5 # half the first-kWh marginal value
elbow_idx = next((i for i, m in enumerate(lp_marg) if m < threshold), len(lp_marg))
if 0 < elbow_idx < len(rows):
elbow = rows[elbow_idx]
print(f"\nLP marginal value drops below 50% of first-kWh value at "
f"capacity ≈ {elbow['capacity_kwh']} kWh "
f"(€{elbow['marginal_lp_eur_per_kwh']:.2f}/kWh added).")
if __name__ == "__main__":
main()