diff --git a/data/raw/extra_batteries.json b/data/raw/extra_batteries.json index fdaecb7..45f1c0b 100644 --- a/data/raw/extra_batteries.json +++ b/data/raw/extra_batteries.json @@ -13,20 +13,5 @@ "recommended": false, "installation": "plugin", "image_url": "" - }, - { - "id": -1, - "title": "Zendure SolarFlow 2400AC 5.8 kWh", - "url": "https://www.solar-outlet.nl/zendure-solarflow-2400ac-thuisbatterijset-58-kwh.html", - "brand": "Zendure", - "brand_slug": "zendure", - "capacity": 5.8, - "power_continuous": 2400, - "price": 1759, - "price_per_kwh": 303.3, - "review_score": null, - "recommended": false, - "installation": "plugin", - "image_url": "" } ] diff --git a/scripts/battery_leaderboard.py b/scripts/battery_leaderboard.py index df8c5e1..e50145c 100644 --- a/scripts/battery_leaderboard.py +++ b/scripts/battery_leaderboard.py @@ -2,11 +2,10 @@ """Run every battery in the catalog through both calculators and rank by payback. Default scenario matches the dynamic-mode store-test: - 4000 kWh/yr net demand, avg retail €0.25, no saldering, 3% inflation. + 4000 kWh/yr demand, avg retail €0.25, no PV, no saldering, 3% inflation. -The simulator works on the raw P1 net-meter trace — whatever PV is on the -roof is already netted in. To model a different PV setup, swap the input -data, not a CLI flag. Output: data/processed/battery_leaderboard.csv. +Override with CLI flags to test a different scenario (e.g., dad's situation +with --pv-kwp 3.0). Output: data/processed/battery_leaderboard.csv. """ from __future__ import annotations @@ -20,6 +19,7 @@ from pluginbattery.sim import ( load_hourly, oracle_daily_schedule, simulate, + synthesize_pv, ) from pluginbattery.store_calc import ( Scenario, @@ -36,6 +36,8 @@ def main() -> None: p.add_argument("--retail", type=float, default=0.25) p.add_argument("--demand", type=float, default=4000.0) p.add_argument("--eta", type=float, default=0.88) + p.add_argument("--pv-kwp", type=float, default=0.0) + p.add_argument("--pv-yield", type=float, default=875.0) p.add_argument("--saldering", action="store_true") p.add_argument("--inflation", type=float, default=0.03) p.add_argument("--top", type=int, default=20, help="How many rows to print") @@ -44,14 +46,15 @@ def main() -> None: # ─── Load and prep the simulation data once ──────────────────────── base = load_hourly("data/raw"); base = apply_nl_tariff(base) df = base.copy() - annual_net_demand = base["raw_demand_kw"].sum() * (8766.0 / len(base)) - df["raw_demand_kw"] = base["raw_demand_kw"] * (args.demand / annual_net_demand) + annual_demand = base["demand_kwh"].sum() * (8766.0 / len(base)) + 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 + if args.pv_kwp > 0: + df = synthesize_pv(df, kwp=args.pv_kwp, + target_kwh_per_kwp_per_year=args.pv_yield) df["export_eur_per_kwh"] = df["eur_per_kwh"] if args.saldering else 0.0 - # Existing-PV export visible at the meter (any hour where raw_demand_kw < 0). - has_existing_pv = bool((df["raw_demand_kw"] < 0).any()) avg_epex = float(df["epex_eur_per_kwh"].mean()) store_params = StoreParams(avg_epex_eur_per_kwh=avg_epex) @@ -77,9 +80,8 @@ def main() -> None: # Cache LP results by (capacity, power) — many batteries share identical specs. lp_cache: dict[tuple[float, float], float] = {} - print(f"Scenario: net demand {args.demand:.0f} kWh/yr, avg retail €{args.retail}, " - f"existing-PV exports={'yes' if has_existing_pv else 'no'}, " - f"saldering={args.saldering}, inflation {args.inflation*100:.1f}%/yr") + print(f"Scenario: demand {args.demand:.0f} kWh/yr, avg retail €{args.retail}, " + f"PV {args.pv_kwp} kWp, saldering={args.saldering}, inflation {args.inflation*100:.1f}%/yr") print(f"Running {len(rows_in)} batteries (caching by (capacity, power))...") print() @@ -102,7 +104,7 @@ def main() -> None: scn = Scenario( capacity_kwh=r["capacity_kwh"], max_charge_kw=r["power_kw"], battery_cost_eur=r["price_eur"], avg_retail_eur_per_kwh=args.retail, - has_pv=has_existing_pv, has_saldering=args.saldering, dynamic_rate=True, + has_pv=args.pv_kwp > 0, has_saldering=args.saldering, dynamic_rate=True, ) sq = store_quote(scn, store_params, inflation=args.inflation) diff --git a/scripts/capacity_sweep.py b/scripts/capacity_sweep.py index 40a6fd0..0f96b7c 100644 --- a/scripts/capacity_sweep.py +++ b/scripts/capacity_sweep.py @@ -48,9 +48,9 @@ def main() -> None: 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)) + annual_demand = base["demand_kwh"].sum() * (8766.0 / len(base)) df = base.copy() - df["raw_demand_kw"] = base["raw_demand_kw"] * (args.demand / annual_demand) + 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 diff --git a/scripts/compare_batteries.py b/scripts/compare_batteries.py index ba5619a..794ec17 100644 --- a/scripts/compare_batteries.py +++ b/scripts/compare_batteries.py @@ -1,10 +1,9 @@ #!/usr/bin/env python3 -"""Run the price-arbitrage oracle for several battery configs and print a comparison. +"""Run the 24h-foresight oracle for several battery configs and print a comparison. Saves per-hour CSVs for each config plus a summary table. Tariff: NL consumer. -The simulator works on the raw P1 net-meter trace — whatever PV is on the -roof is already netted in. Existing-PV exports show up as negative -``raw_demand_kw`` and are absorbed by the battery before crossing the meter. +With --pv-kwp > 0, synthesizes PV output from horizontal irradiance and treats +saldering as full export credit at consumer price (NL pre-2027). """ from __future__ import annotations @@ -19,6 +18,7 @@ from pluginbattery.sim import ( load_hourly, oracle_daily_schedule, simulate, + synthesize_pv, ) CONFIGS = [ @@ -79,6 +79,14 @@ def payback_years(year1_savings: float, cost: float, inflation: float = 0.0) -> def main() -> None: parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--pv-kwp", type=float, default=0.0, + help="PV system size in kWp. 0 = no solar (default).", + ) + parser.add_argument( + "--pv-target", type=float, default=900.0, + help="Calibrate PV to this annual yield in kWh/kWp/year. Default 900 = NL norm.", + ) parser.add_argument( "--saldering", choices=["none", "full"], default="none", help="'full' = export at consumer price (NL pre-2027). " @@ -113,12 +121,11 @@ def main() -> None: elif args.price_mult != 1.0: df["eur_per_kwh"] = df["eur_per_kwh"] * args.price_mult - # Existing PV (if any) is already in raw_demand_kw via the meter signal. - annual_export = float(-df["raw_demand_kw"].clip(upper=0).sum() * 8766.0 / len(df)) - if annual_export > 1.0: - pv_note = f"{annual_export:.0f} kWh/yr exported at meter (existing PV)" - else: - pv_note = "no measurable existing PV in the data" + pv_note = "no solar" + if args.pv_kwp > 0: + df = synthesize_pv(df, kwp=args.pv_kwp, target_kwh_per_kwp_per_year=args.pv_target) + pv_kwh_year = float(df["pv_kwh"].sum() * 8766.0 / len(df)) + pv_note = f"{args.pv_kwp:.1f} kWp PV, calibrated to {pv_kwh_year:.0f} kWh/yr" if args.export_rate is not None: df["export_eur_per_kwh"] = float(args.export_rate) @@ -139,10 +146,10 @@ def main() -> None: print(f"Tariff: import = {price_note}; {sald_note}") if args.price_inflation > 0: print(f"Payback uses {args.price_inflation*100:.1f}%/yr energy-price inflation") - print(f"Meter : {pv_note}\n") + print(f"Solar : {pv_note}\n") # Headline figures with no battery, just for context. - g_no = df["raw_demand_kw"] + g_no = df["demand_kwh"] - df.get("pv_kwh", 0) imp_p = df["eur_per_kwh"] exp_p = df.get("export_eur_per_kwh", df["eur_per_kwh"]) cost_no_battery_total = float((g_no.where(g_no > 0, 0) * imp_p + @@ -180,6 +187,7 @@ def main() -> None: "max_charge_kw": battery.max_charge_kw, "max_discharge_kw": battery.max_discharge_kw, "price_eur": price, + "pv_kwp": args.pv_kwp, "saldering": args.saldering, "annual_savings_eur": round(savings, 2), "kwh_shifted": round(kwh_shifted, 1), diff --git a/scripts/compare_with_store.py b/scripts/compare_with_store.py index 4d793ea..7634fa8 100644 --- a/scripts/compare_with_store.py +++ b/scripts/compare_with_store.py @@ -10,8 +10,8 @@ Examples: python scripts/compare_with_store.py --capacity 5.12 --power 0.8 --cost 1339 --retail 0.25 python scripts/compare_with_store.py --capacity 5.12 --power 2.5 --cost 1339 --retail 0.25 - # Fixed retail, no saldering: - python scripts/compare_with_store.py --fixed --retail 0.28 + # PV scenario, fixed retail, no saldering: + python scripts/compare_with_store.py --pv-kwp 3.0 --fixed --retail 0.28 """ from __future__ import annotations @@ -23,6 +23,7 @@ from pluginbattery.sim import ( load_hourly, oracle_daily_schedule, simulate, + synthesize_pv, ) from pluginbattery.store_calc import ( Scenario, @@ -46,6 +47,10 @@ def main() -> None: help="Fixed-rate mode (no time-of-day variation; kills arbitrage)") p.add_argument("--saldering", action="store_true", help="Full saldering on (export = retail). Default: no saldering.") + # Solar + p.add_argument("--pv-kwp", type=float, default=0.0, help="PV size in kWp; 0 = no PV") + p.add_argument("--pv-yield", type=float, default=875.0, + help="Annual PV yield in kWh/kWp (store uses 875 in their footer)") # Inflation p.add_argument("--inflation", type=float, default=0.03) args = p.parse_args() @@ -54,10 +59,9 @@ def main() -> None: base = load_hourly("data/raw"); base = apply_nl_tariff(base) df = base.copy() - # Demand scaling — scale the entire signed signal so existing-PV exports - # scale with the household. - annual_net_demand = base["raw_demand_kw"].sum() * (8766.0 / len(base)) - df["raw_demand_kw"] = base["raw_demand_kw"] * (args.demand / annual_net_demand) + # Demand scaling. + annual_demand = base["demand_kwh"].sum() * (8766.0 / len(base)) + df["demand_kwh"] = base["demand_kwh"] * (args.demand / annual_demand) # Retail price. if args.fixed: @@ -69,14 +73,16 @@ def main() -> None: df["epex_eur_per_kwh"] = base["epex_eur_per_kwh"] * scale epex_in_use = df["epex_eur_per_kwh"] + # PV. + if args.pv_kwp > 0: + df = synthesize_pv(df, kwp=args.pv_kwp, target_kwh_per_kwp_per_year=args.pv_yield) + # Saldering. if args.saldering: df["export_eur_per_kwh"] = df["eur_per_kwh"] else: df["export_eur_per_kwh"] = 0.0 # post-saldering default in this script - has_existing_pv = bool((df["raw_demand_kw"] < 0).any()) - bat = Battery( capacity_kwh=args.capacity, max_charge_kw=args.power, @@ -96,7 +102,7 @@ def main() -> None: max_charge_kw=args.power, battery_cost_eur=args.cost, avg_retail_eur_per_kwh=args.retail, - has_pv=has_existing_pv, + has_pv=args.pv_kwp > 0, has_saldering=args.saldering, dynamic_rate=not args.fixed, ) @@ -108,8 +114,9 @@ def main() -> None: print(f" Tariff : avg retail €{args.retail}, " f"{'FIXED' if args.fixed else 'DYNAMIC'} rate, " f"{'WITH' if args.saldering else 'NO'} saldering") - print(f" PV : {'existing PV in meter signal' if has_existing_pv else 'no measurable existing PV'}") - print(f" Net demand : {args.demand:.0f} kWh/yr") + print(f" PV : {args.pv_kwp} kWp" + f"{' @ ' + str(args.pv_yield) + ' kWh/kWp/yr' if args.pv_kwp else ''}") + print(f" Demand : {args.demand:.0f} kWh/yr") print(f" Inflation : {args.inflation*100:.1f}%/yr") print() print(f"{'':22s} {'year-1 €':>10s} {'payback':>9s}") diff --git a/scripts/import_dad_p1.py b/scripts/import_dad_p1.py index 2c114e2..2595478 100644 --- a/scripts/import_dad_p1.py +++ b/scripts/import_dad_p1.py @@ -6,9 +6,9 @@ meter readings, every 15 min, 2025 calendar year, Europe/Amsterdam local time). Output: data/raw/p1_hourly.csv (timestamp UTC, power_w mean over hour). Net of PV exports — `power_w` can be negative when his solar exceeded household -demand. The simulator interprets `raw_demand_kw = power_w / 1000`, and a -plug-in battery's discharge cap `max(0, raw_demand_kw)` correctly forbids -discharge during net-export hours (the meter is already flowing the wrong way). +demand. The simulator interprets `demand_kwh = power_w / 1000`, and a plug-in +battery's discharge cap `max(0, demand - pv)` correctly handles negative +demand by forbidding discharge during surplus hours. The simulator's window is 2023-09-01 → 2024-09-01 UTC. Dad's data is the 2025 calendar year. We align by (month, day, hour) so seasonal patterns diff --git a/scripts/timeline_inspect.py b/scripts/timeline_inspect.py index a7c6838..0cdc337 100644 --- a/scripts/timeline_inspect.py +++ b/scripts/timeline_inspect.py @@ -46,8 +46,8 @@ def main() -> None: # ─── Build the same df the web app sees ─── base = apply_nl_tariff(load_hourly("data/raw")) df = base.copy() - annual = base["raw_demand_kw"].sum() * 8766 / len(base) - df["raw_demand_kw"] = base["raw_demand_kw"] * (args.demand_kwh / annual) + annual = base["demand_kwh"].sum() * 8766 / len(base) + df["demand_kwh"] = base["demand_kwh"] * (args.demand_kwh / annual) 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 @@ -62,10 +62,10 @@ def main() -> None: out = out.assign( import_kwh=np.maximum(0.0, g), export_kwh=np.maximum(0.0, -g), - no_bat_grid=out["raw_demand_kw"], # equals the original net since pv=0 - cost_no_battery_eur=np.where(out["raw_demand_kw"] > 0, - out["raw_demand_kw"] * out["eur_per_kwh"], - out["raw_demand_kw"] * out["export_eur_per_kwh"]), + no_bat_grid=out["demand_kwh"], # equals the original net since pv=0 + cost_no_battery_eur=np.where(out["demand_kwh"] > 0, + out["demand_kwh"] * out["eur_per_kwh"], + out["demand_kwh"] * out["export_eur_per_kwh"]), cost_with_battery_eur=np.where(g > 0, g * out["eur_per_kwh"], g * out["export_eur_per_kwh"]), ) @@ -89,7 +89,7 @@ def main() -> None: print(f" {'hr':>2s} {'price':>5s} {'no-bat meter':>15s} " f"{'charge':>7s} {'disch':>7s} {'SoC':>5s} {'with-bat meter':>16s} {'sav':>6s}") for ts, r in day.iterrows(): - net_no = r["raw_demand_kw"] + net_no = r["demand_kwh"] mark = "↓ exp" if net_no < -0.01 else "↑ imp" if net_no > 0.01 else "·" grid = r["grid_kwh_with_battery"] grid_mark = "↓ exp" if grid < -0.01 else "↑ imp" if grid > 0.01 else "·" @@ -97,8 +97,8 @@ def main() -> None: f"{net_no:+7.2f} {mark:5s} " f"{r['charge_kwh']:>6.2f} {r['discharge_kwh']:>6.2f} " f"{r['soc_kwh']:>4.2f} {grid:+8.2f} {grid_mark:5s} €{r['savings_eur']:+5.2f}") - d_imp_no = max(0.0, day['raw_demand_kw'].clip(lower=0).sum()) - d_exp_no = -day['raw_demand_kw'].clip(upper=0).sum() + d_imp_no = max(0.0, day['demand_kwh'].clip(lower=0).sum()) + d_exp_no = -day['demand_kwh'].clip(upper=0).sum() d_imp_yes = day['import_kwh'].sum() d_exp_yes = day['export_kwh'].sum() d_sav = day['savings_eur'].sum() @@ -113,8 +113,8 @@ def main() -> None: week_savings = win["savings_eur"].sum() week_charge = win["charge_kwh"].sum() week_disch = win["discharge_kwh"].sum() - week_surplus = (-win["raw_demand_kw"].clip(upper=0)).sum() - week_imports = win["raw_demand_kw"].clip(lower=0).sum() + week_surplus = (-win["demand_kwh"].clip(upper=0)).sum() + week_imports = win["demand_kwh"].clip(lower=0).sum() annual_savings = out["savings_eur"].sum() annualised = week_savings * 365.25 / args.days diff --git a/src/pluginbattery/sim.py b/src/pluginbattery/sim.py index 02fa178..b8e3471 100644 --- a/src/pluginbattery/sim.py +++ b/src/pluginbattery/sim.py @@ -1,22 +1,13 @@ """Hourly home-battery simulator. -The simulator operates on a single signed signal: ``raw_demand_kw``, the P1 -net-meter reading. Positive = importing, negative = exporting (e.g. from -existing rooftop PV already netted at the meter). We do **not** model solar -as a separate signal — whatever the meter shows is the input. To ask "what -if the PV were different" you need actual measurements of that situation; -synthesizing PV from a remote irradiance sensor was tried and didn't survive -contact with reality (different roof, different orientation, different cloud -cover). - Three pieces: - load_hourly(): CSVs in -> single hourly DataFrame, missing-hour rows dropped. - simulate(df, bat, sched): walks SoC hour by hour, clamps to physical limits, returns trajectory + grid flows + costs. - - oracle_daily_schedule(): greedy self-consumption (Pass 1) + perfect-foresight - per-day price arbitrage (Pass 2). + - oracle_daily_schedule(): solves a 24h LP per UTC day for perfect-foresight + within-day arbitrage. No cross-day shifting. -Anything that decides what to do returns a (n, 2) schedule of +Anything that decides what to do (oracle, rules, ML) returns a (n, 2) schedule of [charge_kwh, discharge_kwh]; simulate() is the single source of truth for state. """ from __future__ import annotations @@ -26,6 +17,7 @@ from pathlib import Path import numpy as np import pandas as pd +from scipy.optimize import linprog @dataclass(frozen=True) @@ -48,6 +40,39 @@ ECOFLOW_STREAM_AC = Battery( ) +def synthesize_pv( + df: pd.DataFrame, + kwp: float, + target_kwh_per_kwp_per_year: float = 900.0, + irradiance_cap_w_m2: float = 1000.0, +) -> pd.DataFrame: + """Synthesize hourly PV output (kWh) from horizontal irradiance. + + Output shape follows the irradiance signal (so timing tracks real local + weather). The total is calibrated to `target_kwh_per_kwp_per_year` — + 900 kWh/kWp/year is the NL norm for a south-facing 35° roof. This + sidesteps the question "what's the right system factor?" and absorbs any + sensor bias automatically: if the irradiance reads 22% low, the implied + factor goes up to compensate. + + Outlier handling: irradiance clipped to [0, irradiance_cap_w_m2]. + """ + out = df.copy() + irr = out["irradiance_w_m2"].clip(lower=0.0, upper=irradiance_cap_w_m2) + + # Annualised horizontal irradiance density for *this* dataset, kWh/m²/year. + hours = len(df) + annual_irradiance_kwh_m2 = (irr.sum() / 1000.0) * (8766.0 / hours) + if annual_irradiance_kwh_m2 <= 0: + raise ValueError("Irradiance data is all zero or negative.") + + system_factor = target_kwh_per_kwp_per_year / annual_irradiance_kwh_m2 + out["pv_kwh"] = irr / 1000.0 * kwp * system_factor + out.attrs["pv_system_factor"] = system_factor + out.attrs["pv_annual_kwh"] = float(out["pv_kwh"].sum() * 8766.0 / hours) + return out + + def apply_nl_tariff( df: pd.DataFrame, vat: float = 0.21, @@ -70,16 +95,17 @@ def apply_nl_tariff( def load_hourly(data_dir: Path | str = "data/raw") -> pd.DataFrame: - """Read prices + P1 CSVs into one hourly DataFrame indexed by UTC timestamp. + """Read prices/P1/solar CSVs into one hourly DataFrame indexed by UTC timestamp. - Drops every hour where either signal is missing (inner join + dropna). - Adds ``raw_demand_kw = power_w / 1000`` — the signed P1 net meter - reading. Negative = exporting, positive = importing. + Drops every hour where any of the three signals is missing (inner join + + dropna). Adds a derived `demand_kwh` column = power_w / 1000 (since power_w + is the hour mean and 1 W * 1 h = 1 Wh). """ data_dir = Path(data_dir) files = { "eur_per_kwh": "prices_hourly.csv", - "power_w": "p1_hourly.csv", + "power_w": "p1_hourly.csv", + "irradiance_w_m2": "solar_hourly.csv", } parts = [] for col, fname in files.items(): @@ -91,7 +117,7 @@ def load_hourly(data_dir: Path | str = "data/raw") -> pd.DataFrame: sub.index = sub.index.tz_convert("UTC") parts.append(sub[[col]]) df = pd.concat(parts, axis=1, join="inner").dropna().sort_index() - df["raw_demand_kw"] = df["power_w"] / 1000.0 + df["demand_kwh"] = df["power_w"] / 1000.0 return df @@ -103,13 +129,13 @@ def simulate(df: pd.DataFrame, battery: Battery, schedule: np.ndarray) -> pd.Dat authoritative state engine, the schedule is only a *request*. Column conventions read by this function: - - raw_demand_kw: signed meter reading (always required). +import, −export. - - eur_per_kwh: import price (always required). + - eur_per_kwh: import price (always required) + - pv_kwh: optional, hourly PV production - export_eur_per_kwh: optional, export credit price. Defaults to eur_per_kwh (full saldering) if missing. - Per hour: ``g = raw_demand_kw + charge − discharge``. Cost uses the - import price when g > 0, the export price when g < 0. + Per hour: `g = demand − pv + charge − discharge`. Cost uses the import + price when g > 0, the export price when g < 0. """ n = len(df) if schedule.shape != (n, 2): @@ -125,7 +151,8 @@ def simulate(df: pd.DataFrame, battery: Battery, schedule: np.ndarray) -> pd.Dat d_arr = np.empty(n) soc_arr = np.empty(n) - raw = df["raw_demand_kw"].to_numpy() + demand = df["demand_kwh"].to_numpy() + pv = df["pv_kwh"].to_numpy() if "pv_kwh" in df.columns else np.zeros(n) requested_c = schedule[:, 0] requested_d = schedule[:, 1] @@ -139,8 +166,8 @@ def simulate(df: pd.DataFrame, battery: Battery, schedule: np.ndarray) -> pd.Dat d = max(0.0, min(requested_d[t], pd_max, avail_kwh_grid)) if not battery.allows_export: # Plug-in: cannot push current backwards through the meter. - # Discharge can't take the net flow below zero. - d = min(d, max(0.0, raw[t])) + # Discharge can't exceed the *net* demand after solar. + d = min(d, max(0.0, demand[t] - pv[t])) soc = soc + c * eta - d / eta # Numerical guard rail for floating point drift. @@ -151,6 +178,8 @@ def simulate(df: pd.DataFrame, battery: Battery, schedule: np.ndarray) -> pd.Dat soc_arr[t] = soc out = df.copy() + if "pv_kwh" not in out.columns: + out["pv_kwh"] = 0.0 out["charge_kwh"] = c_arr out["discharge_kwh"] = d_arr out["soc_kwh"] = soc_arr @@ -161,7 +190,7 @@ def simulate(df: pd.DataFrame, battery: Battery, schedule: np.ndarray) -> pd.Dat else: export_price = import_price # default = full saldering - g_no = raw + g_no = (out["demand_kwh"] - out["pv_kwh"]).to_numpy() g_yes = g_no + c_arr - d_arr out["grid_kwh_no_battery"] = g_no @@ -172,12 +201,32 @@ def simulate(df: pd.DataFrame, battery: Battery, schedule: np.ndarray) -> pd.Dat return out -def _pass1_reactive(df: pd.DataFrame, battery: Battery) -> tuple[np.ndarray, np.ndarray]: - """Greedy self-consumption: charge from surplus, discharge to cover demand. +def oracle_daily_schedule(df: pd.DataFrame, battery: Battery) -> np.ndarray: + """Two-pass dispatch: greedy self-consumption + daily grid arbitrage. - Pure function of raw_demand_kw. No foresight — exactly the "default - mode" behaviour of every plug-in battery firmware out there. - Returns (schedule, soc) where soc has length n+1 (soc[0] = initial). + This matches what a modern plug-in battery with smart-charging firmware + (Tibber controllers, Marstek 'Smart Charging' mode, EcoFlow with dynamic + tariff integration) actually does in practice — and it always produces + visually intuitive schedules without LP indifference artefacts. + + Pass 1 (greedy self-consumption): + For each hour in order, charge from any meter export, discharge to + cover any meter import. Exactly the "default mode" behaviour. + + Pass 2 (per-UTC-day arbitrage on residual capacity): + For each day, repeatedly find the most profitable + cheap-hour → expensive-hour pair, considering: + - Both charge and discharge already have residual room (power + and SoC headroom). + - Discharge fits within demand at the expensive hour (plug-in + cannot push to grid; allows_export=True relaxes this). + - Profit per shifted kWh = η · price[discharge] − price[charge] / η + must be positive. + Execute the most profitable, update SoC trajectory, repeat. Stops + when no profitable cycle remains for that day. + + Pass 1 is order-of-day dependent (real firmware acts moment-by-moment), + so it always runs first. Pass 2 sees what's left. """ cap = battery.capacity_kwh pc_max = battery.max_charge_kw @@ -185,50 +234,38 @@ def _pass1_reactive(df: pd.DataFrame, battery: Battery) -> tuple[np.ndarray, np. eta = battery.round_trip_eff ** 0.5 n = len(df) - raw = df["raw_demand_kw"].to_numpy() + demand = df["demand_kwh"].to_numpy() + pv = df["pv_kwh"].to_numpy() if "pv_kwh" in df.columns else np.zeros(n) + net = demand - pv + prices = df["eur_per_kwh"].to_numpy() schedule = np.zeros((n, 2)) soc = np.zeros(n + 1) soc[0] = battery.initial_soc_kwh + # ── Pass 1: greedy self-consumption ───────────────────────────────── for t in range(n): - if raw[t] < -1e-9: # exporting (surplus from existing PV) - surplus = -raw[t] + if net[t] < -1e-9: # exporting + surplus = -net[t] room = (cap - soc[t]) / eta if eta > 0 else 0.0 c = max(0.0, min(pc_max, surplus, room)) schedule[t, 0] = c soc[t + 1] = min(cap, soc[t] + c * eta) - elif raw[t] > 1e-9: # importing + elif net[t] > 1e-9: # importing avail = soc[t] * eta - cap_d = pd_max if battery.allows_export else raw[t] + cap_d = pd_max if battery.allows_export else net[t] d = max(0.0, min(pd_max, avail, cap_d)) schedule[t, 1] = d soc[t + 1] = max(0.0, soc[t] - d / eta) else: soc[t + 1] = soc[t] - return schedule, soc - -def _pass2_arbitrage(df: pd.DataFrame, battery: Battery, - schedule: np.ndarray, soc: np.ndarray) -> None: - """Per-UTC-day arbitrage layered on top of (schedule, soc). Mutates in place. - - For each day, repeatedly find the most profitable (ci, di) pair where - profit = η·price[di] − price[ci]/η > 0 and both endpoints have power + - SoC headroom. Execute, update soc trajectory, repeat. Stops when no - profitable cycle remains. - """ - cap = battery.capacity_kwh - pc_max = battery.max_charge_kw - pd_max = battery.max_discharge_kw - eta = battery.round_trip_eff ** 0.5 + # ── Pass 2: per-day arbitrage on residual capacity ────────────────── if eta <= 0: - return + return schedule inv_eta = 1.0 / eta - raw = df["raw_demand_kw"].to_numpy() - prices = df["eur_per_kwh"].to_numpy() - + # Index hours by UTC date (timeline is short enough that this is cheap). by_date: dict = {} for i, ts in enumerate(df.index): by_date.setdefault(ts.date(), []).append(i) @@ -242,6 +279,8 @@ def _pass2_arbitrage(df: pd.DataFrame, battery: Battery, room_c = pc_max - schedule[ci, 0] if room_c < 1e-9: continue + # SoC headroom across (ci, end-of-day]: how much can we add + # to soc[ci+1..last+1] without exceeding cap? tail = soc[ci + 1 : last + 2] hr = (cap - tail.max()) / eta if tail.size else 0.0 if hr < 1e-9: @@ -252,16 +291,25 @@ def _pass2_arbitrage(df: pd.DataFrame, battery: Battery, if room_d < 1e-9: continue if not battery.allows_export: - # Plug-in: extra discharge limited by the post-Pass-1 import demand. - room_d = min(room_d, raw[di] - schedule[di, 1]) + # Discharge cannot exceed demand at di. + room_d = min(room_d, net[di] - schedule[di, 1]) if room_d < 1e-9: continue + # SoC must remain ≥ 0 across (ci, di]; new schedule has + # +c at ci shifting trajectory up, then -d at di pulling + # it down. Trajectory between ci+1..di stays elevated by + # +x·η; minimum allowed kWh limited by min SoC in that + # range plus future-discharge availability. profit = eta * prices[di] - prices[ci] * inv_eta if profit <= 1e-6: continue + + # Max charge in AC kWh: limited by power, capacity headroom, + # discharge-side capacity. max_kwh = min(room_c, hr, room_d * inv_eta) if max_kwh < 1e-6: continue + if best is None or profit > best[2]: best = (ci, di, profit, max_kwh) if best is None: @@ -270,20 +318,130 @@ def _pass2_arbitrage(df: pd.DataFrame, battery: Battery, ci, di, _, kwh = best schedule[ci, 0] += kwh schedule[di, 1] += kwh * eta + # Trajectory: SoC between ci+1 and di stays elevated by kwh*eta. + # SoC at di+1 onwards returns to original (+kwh*eta) - (kwh*eta) = unchanged). for t in range(ci, di): soc[t + 1] += kwh * eta - # soc[di+1] onwards: +kwh*eta from charge cancels with -kwh*eta - # from discharge → unchanged. + # Beyond di, soc[di+1] = previous - d/eta = previous - kwh + # We added kwh*eta at soc[di], then removed kwh*eta on the + # discharge so soc[di+1] should already match — but we updated + # the loop with soc trajectory; need to ensure di+1 onward stays + # consistent. + # Easier: walk soc again from di onward to verify consistency. + # Since both sides cancel, soc[di+1..end] should be unchanged. + + return schedule + + +def _oracle_daily_schedule_legacy(df: pd.DataFrame, battery: Battery) -> np.ndarray: + """Legacy per-day LP (daily-myopic). Kept for reference / regression + testing only — the public oracle_daily_schedule() now solves the whole + window in one go to avoid SoC-resets-to-zero-overnight artefacts.""" + cap = battery.capacity_kwh + pc_max = battery.max_charge_kw + pd_max = battery.max_discharge_kw + eta = battery.round_trip_eff ** 0.5 + + n = len(df) + schedule = np.zeros((n, 2)) + soc = battery.initial_soc_kwh + + # Position lookup so we can fill the global schedule from each day's solution. + pos = {ts: i for i, ts in enumerate(df.index)} + + has_pv = "pv_kwh" in df.columns + has_export_price = "export_eur_per_kwh" in df.columns + + for _, day_df in df.groupby(df.index.date, sort=True): + m = len(day_df) + import_price = day_df["eur_per_kwh"].to_numpy() + export_price = ( + day_df["export_eur_per_kwh"].to_numpy() if has_export_price else import_price + ) + demand = day_df["demand_kwh"].to_numpy() + pv = day_df["pv_kwh"].to_numpy() if has_pv else np.zeros(m) + + # Variables: x = [c_0..c_{m-1}, d_0..d_{m-1}, gI_0..gI_{m-1}, gE_0..gE_{m-1}] + # c = charge from bus to battery (kWh AC) + # d = discharge from battery to bus (kWh AC) + # gI = grid import to bus (kWh, ≥ 0) + # gE = grid export from bus (kWh, ≥ 0) + # Energy balance per hour: gI - gE = demand - pv + c - d + # → gI - gE - c + d = demand - pv + # Cost = sum(import_price * gI - export_price * gE). + # When import_price > export_price, the LP automatically picks + # min(gI, gE) = 0 because keeping both > 0 costs (imp - exp) > 0. + # + # Tiny tie-breakers: when several surplus hours are equally good for + # charging (true whenever export_price is constant and we're capturing + # surplus), the simplex would otherwise pick an arbitrary subset. + # Adding a microscopic penalty proportional to hour-index makes early + # charging strictly cheaper than late charging by a fraction of a + # cent — so the battery visibly fills as soon as surplus appears. + # Symmetric for discharge: prefer late hours (when prices are + # typically highest anyway, but ties happen during evening plateaus). + EPS = 1e-6 + hr = np.arange(m) + f = np.concatenate([ + EPS * hr, # c: prefer earlier charging (smaller h ⇒ less cost) + EPS * (m - 1 - hr), # d: prefer later discharging + import_price, # pay for imports + -export_price, # earn export credit + ]) + + # Variable bounds. + c_bounds = [(0.0, pc_max)] * m + if battery.allows_export: + d_bounds = [(0.0, pd_max)] * m + else: + # Plug-in: discharge ≤ net demand after solar (battery cannot push current backwards). + d_bounds = [ + (0.0, min(pd_max, max(0.0, demand[t] - pv[t]))) for t in range(m) + ] + # Physical caps on grid flow per hour. Without these the LP can find + # zero-cost cycles when export_price ≈ 0 (charge from cheap surplus, + # discharge into cheap export — both sides zero, so any volume is "free"). + # gI: practical max draw is demand + max charge; cap generously at 100 kWh/h. + # gE: bounded by *natural* export capacity = max(0, pv − demand). This + # also handles the case where the demand series is already net of an + # external PV (so demand can be negative even when pv = 0). + gi_bounds = [(0.0, 100.0)] * m + natural_export = np.maximum(0.0, pv - demand) + if battery.allows_export: + ge_bounds = [(0.0, max(1e-9, natural_export[t] + pd_max)) for t in range(m)] + else: + ge_bounds = [(0.0, max(1e-9, natural_export[t])) for t in range(m)] + bounds = c_bounds + d_bounds + gi_bounds + ge_bounds + + # Energy balance equality (m rows): -c[t] + d[t] + gI[t] - gE[t] = demand[t] - pv[t] + I = np.eye(m) + A_eq = np.hstack([-I, I, I, -I]) + b_eq = demand - pv + + # SoC inequality (2m rows). SoC(t) = soc0 + eta*cumsum(c)[t] - cumsum(d)[t]/eta. + # Apply only to c, d columns; gI, gE are zero in these rows. + L = np.tril(np.ones((m, m))) + Z = np.zeros((m, m)) + A_ub = np.vstack([ + np.hstack([-eta * L, (1.0 / eta) * L, Z, Z]), # -SoC ≤ soc0 ⇒ SoC ≥ 0 + np.hstack([ eta * L, -(1.0 / eta) * L, Z, Z]), # SoC ≤ cap - soc0 + soc0 + ]) + b_ub = np.concatenate([np.full(m, soc), np.full(m, cap - soc)]) + + result = linprog( + f, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="highs" + ) + if not result.success: + raise RuntimeError(f"LP failed for {day_df.index[0].date()}: {result.message}") + + c_day = result.x[:m] + d_day = result.x[m : 2 * m] + + for t, ts in enumerate(day_df.index): + schedule[pos[ts], 0] = c_day[t] + schedule[pos[ts], 1] = d_day[t] + + soc = soc + c_day.sum() * eta - d_day.sum() / eta + soc = max(0.0, min(cap, soc)) - -def oracle_daily_schedule(df: pd.DataFrame, battery: Battery) -> np.ndarray: - """Greedy self-consumption (Pass 1) + perfect-foresight per-day arbitrage (Pass 2). - - The ceiling for what a controller could achieve given a day-ahead price - series. Cross-day shifting is intentionally not allowed (the next day's - SoC starts wherever the previous day left it, but no within-day - arbitrage looks across the midnight boundary). - """ - schedule, soc = _pass1_reactive(df, battery) - _pass2_arbitrage(df, battery, schedule, soc) return schedule diff --git a/src/pluginbattery/static/app.js b/src/pluginbattery/static/app.js index bb598f5..6a1f4fe 100644 --- a/src/pluginbattery/static/app.js +++ b/src/pluginbattery/static/app.js @@ -56,6 +56,8 @@ function readScenario() { return { demand_kwh: parseFloat(fd.get("demand_kwh")), retail: parseFloat(fd.get("retail")), + pv_kwp: parseFloat(fd.get("pv_kwp")), + pv_yield: 875.0, fixed_rate: fd.get("fixed_rate") === "true", saldering: fd.get("saldering") === "true", eta: parseFloat(fd.get("eta")), @@ -265,7 +267,7 @@ const tlSummary = mustFind("tl-summary"); let activeBattery = null; // the row clicked let tlData = null; // last fetched timeline payload -let charts = { price: null, grid: null, soc: null }; +let charts = { irr: null, soc: null, grid: null }; function openModalForBattery(row) { activeBattery = row; @@ -290,7 +292,7 @@ function openModalForBattery(row) { function closeModal() { modal.hidden = true; for (const c of Object.values(charts)) c?.destroy?.(); - charts = { price: null, grid: null, soc: null }; + charts = { irr: null, soc: null, grid: null }; } modal.querySelectorAll('[data-close="modal"]').forEach( @@ -323,7 +325,7 @@ async function refreshTimeline() { start: tlDate.value || "2024-06-17", days: tlDays.value, demand_kwh: sc.demand_kwh, retail: sc.retail, - fixed_rate: sc.fixed_rate, + pv_kwp: sc.pv_kwp, fixed_rate: sc.fixed_rate, saldering: sc.saldering, eta: sc.eta, }); if (sc.export_rate !== null) params.set("export_rate", sc.export_rate); @@ -354,12 +356,12 @@ function paintTimeline(d) { tlSummary.replaceChildren(); for (const [k, v] of [ ["window cost", `€${cost.toFixed(2)}`], - ["imports", `${imports.toFixed(2)} kWh`], - ["exports", `${exports_.toFixed(2)} kWh`], + ["imports", `${imports.toFixed(1)} kWh`], + ["exports", `${exports_.toFixed(1)} kWh`], ["battery throughput", batteryOn - ? `${s.battery_charge_kwh.toFixed(2)} in / ${s.battery_discharge_kwh.toFixed(2)} out kWh` + ? `${s.battery_charge_kwh.toFixed(1)} in / ${s.battery_discharge_kwh.toFixed(1)} out kWh` : `— (battery off)`], - ["window saving", `€${savings.toFixed(4)}`], + ["window saving", `€${savings.toFixed(2)}`], ]) { const row = document.createElement("div"); const lbl = document.createElement("span"); lbl.textContent = k; @@ -368,26 +370,26 @@ function paintTimeline(d) { } // X axis: unix seconds (uPlot wants epoch seconds) - const xs = d.hours.map(h => Date.parse(h.ts) / 1000); - const price = d.hours.map(h => h.price); - const soc = batteryOn ? d.hours.map(h => h.soc) : d.hours.map(() => 0); - const grid = batteryOn ? d.hours.map(h => h.grid_with_bat) - : d.hours.map(h => h.grid_no_bat); + const xs = d.hours.map(h => Date.parse(h.ts) / 1000); + const irr = d.hours.map(h => h.irradiance); + const soc = batteryOn ? d.hours.map(h => h.soc) : d.hours.map(() => 0); + const grid = batteryOn ? d.hours.map(h => h.grid_with_bat) + : d.hours.map(h => h.grid_no_bat); - drawChart("price", ["Time", "€/kWh"], xs, price, "tl-chart-price", { - stroke: "#a78bfa", fill: "rgba(167, 139, 250, 0.2)", - fmtY: (v) => v == null ? "" : `€${v.toFixed(3)}/kWh`, + drawChart("irr", ["Time", "W/m²"], xs, irr, "tl-chart-irradiance", { + stroke: "#fbbf24", fill: "rgba(251, 191, 36, 0.2)", + fmtY: (v) => v == null ? "" : `${v.toFixed(0)} W/m²`, }); - drawChart("grid", ["Time", "kWh/h"], xs, grid, "tl-chart-grid", { + drawChart("soc", ["Time", "kWh"], xs, soc, "tl-chart-soc", { + stroke: "#4ade80", fill: "rgba(74, 222, 128, 0.2)", + fmtY: (v) => v == null ? "" : `${v.toFixed(2)} kWh`, + }); + drawChart("grid", ["Time", "kWh/h"], xs, grid, "tl-chart-grid", { stroke: batteryOn ? "#93c5fd" : "#fb7185", fill: batteryOn ? "rgba(147, 197, 253, 0.18)" : "rgba(251, 113, 133, 0.18)", zeroline: true, fmtY: (v) => v == null ? "" : `${v >= 0 ? "+" : ""}${v.toFixed(2)} kWh/h`, }); - drawChart("soc", ["Time", "kWh"], xs, soc, "tl-chart-soc", { - stroke: "#4ade80", fill: "rgba(74, 222, 128, 0.2)", - fmtY: (v) => v == null ? "" : `${v.toFixed(2)} kWh`, - }); } function drawChart(key, axes, xs, ys, containerId, opts = {}) { diff --git a/src/pluginbattery/static/style.css b/src/pluginbattery/static/style.css index 0c314f5..f3f4781 100644 --- a/src/pluginbattery/static/style.css +++ b/src/pluginbattery/static/style.css @@ -238,25 +238,16 @@ table#leaderboard tbody tr a { color: inherit; text-decoration: underline dotted } .modal-summary div { display: flex; justify-content: space-between; gap: 0.5rem; } .modal-summary div strong { color: var(--accent); } -.modal-strategy { - border: 1px solid var(--border); border-radius: 8px; background: var(--panel2); - padding: 0.45rem 0.75rem 0.6rem; margin-bottom: 0.75rem; - display: flex; flex-wrap: wrap; gap: 0.25rem 1.25rem; align-items: center; -} -.modal-strategy legend { padding: 0 0.4rem; font-size: 0.78rem; color: var(--muted); } -.modal-strategy .radio { display: inline-flex; align-items: center; gap: 0.4rem; - font-size: 0.85rem; color: var(--text); cursor: pointer; } -.modal-strategy .radio input { margin: 0; } .modal-charts { display: flex; flex-direction: column; gap: 0.6rem; } .modal-charts figure { margin: 0; background: var(--panel); border: 1px solid var(--border); - border-radius: 8px; padding: 0.6rem 0.75rem 1rem; } + border-radius: 8px; padding: 0.6rem 0.75rem 0.4rem; } .modal-charts figcaption { font-size: 0.78rem; color: var(--muted); margin-bottom: 0.3rem; } -.modal-charts .chart { width: 100%; height: 150px; } -.modal-charts figure:nth-child(4) .chart { height: 230px; } /* taller grid panel (4th now) */ +.modal-charts .chart { width: 100%; height: 110px; } +.modal-charts figure:nth-child(3) .chart { height: 200px; } /* taller grid panel */ /* uPlot dark-theme tweaks */ .uplot { color: var(--text); } -.uplot .u-legend { color: var(--text); font-size: 0.75rem; padding: 1rem 0 0.6rem; } +.uplot .u-legend { color: var(--text); font-size: 0.75rem; padding: 0.15rem 0; } .uplot .u-legend th, .uplot .u-legend td { color: var(--text); padding: 0 0.4rem; } .uplot .u-legend .u-marker { display: inline-block; } .uplot .u-axis { color: var(--muted); } diff --git a/src/pluginbattery/templates/index.html b/src/pluginbattery/templates/index.html index 6b3c90a..dcce9cf 100644 --- a/src/pluginbattery/templates/index.html +++ b/src/pluginbattery/templates/index.html @@ -5,7 +5,7 @@ Honest battery payback — vs thuisbatterijgids.nl - + @@ -22,18 +22,21 @@