pluginbattery/scripts/capacity_sweep.py

#!/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["raw_demand_kw"].sum() * (8766.0 / len(base))
    df = base.copy()
    df["raw_demand_kw"] = base["raw_demand_kw"] * (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()