diff --git a/scripts/analysis/battery_sizing_screen.py b/scripts/analysis/battery_sizing_screen.py new file mode 100644 index 0000000..23eb50d --- /dev/null +++ b/scripts/analysis/battery_sizing_screen.py @@ -0,0 +1,368 @@ +#!/usr/bin/env python3 +""" +Ekonomický screening velikosti baterie (15min, jednobusová energie). + +Typicky BA81: fixní nákup + prodej spot, limity výkonu z baterie min(0,5C, střídač), +export/import podle připojení. Načte OTE z Postgres (stejná DB jako EMS) nebo z CSV. + +Příklad: + python3 scripts/analysis/battery_sizing_screen.py \\ + --db \\ + --date-from 2024-04-01 --date-to 2026-04-01 \\ + --battery-kwh 12.5 32 48 \\ + --load-kw 1.2 \\ + --pv-daily-kwh-summer 55 --pv-daily-kwh-winter 12 \\ + --sell-margin-fixed -0.02 \\ + --buy-vat-kwh 4.443 \\ + --capex-per-kwh 9000 + +Vyžaduje: pip install pulp (volitelně psycopg2 pro --db). + +Omezení modelu: syntetický denní tvar FVE (kalibruj --pv-daily-kwh-* podle měření); +mikroinvertory / GEN nejsou; zelený bonus není v účelové funkci; nákup je jedna flat +sazba vč. DPH (reálné NT/VT přes HDO přidej později). Výsledek = screening, ne nabídka. +""" +from __future__ import annotations + +import argparse +import csv +import math +import os +import sys +from dataclasses import dataclass +from datetime import date, datetime, timedelta +from typing import Iterable, Sequence + +try: + import pulp +except ImportError: + print("Instaluj PuLP: pip install pulp", file=sys.stderr) + raise + +DT_H = 0.25 # 15 min +SLOTS_PER_DAY = 96 + + +@dataclass +class SiteLimits: + max_export_w: float = 16_000.0 + max_import_w: float = 17_000.0 + inv_batt_max_w: float = 12_000.0 # strop střídače z baterie / nabíjení + c_rate: float = 0.5 # P_batt = min(c_rate * E_kWh * 1000, inv_batt_max_w) + eta_charge: float = 0.95 + eta_discharge: float = 0.95 + soc_min_frac: float = 0.10 + soc_max_frac: float = 0.95 + + +def batt_power_cap_w(usable_kwh: float, site: SiteLimits) -> float: + return min(site.c_rate * usable_kwh * 1000.0, site.inv_batt_max_w) + + +def summer_day(d: date) -> bool: + m = d.month + return m >= 4 and m <= 9 + + +def pv_shape_96() -> list[float]: + """Nenormalizovaný denní tvar (96 slotů), plocha = 1 po normalizaci.""" + w = [0.0] * SLOTS_PER_DAY + for t in range(SLOTS_PER_DAY): + h = t / 4.0 # hodiny od půlnoci + if 5.5 <= h <= 20.5: + w[t] = max(0.0, math.sin(math.pi * (h - 5.5) / 15.0)) ** 1.2 + else: + w[t] = 0.0 + s = sum(w) + if s <= 0: + return [1.0 / SLOTS_PER_DAY] * SLOTS_PER_DAY + return [x / s for x in w] + + +def daily_pv_wh(d: date, summer_kwh: float, winter_kwh: float, shape: Sequence[float]) -> list[float]: + base = summer_kwh if summer_day(d) else winter_kwh + return [base * 1000.0 * sh for sh in shape] + + +def daily_load_wh(load_kw: float) -> list[float]: + e_per_slot = load_kw * 1000.0 * DT_H + return [e_per_slot] * SLOTS_PER_DAY + + +def effective_sell_kc_kwh(raw_ote: float, margin_fixed: float, margin_pct: float) -> float: + return raw_ote + margin_fixed + (raw_ote * margin_pct / 100.0) + + +def load_prices_csv(path: str) -> list[tuple[datetime, float]]: + out: list[tuple[datetime, float]] = [] + with open(path, newline="", encoding="utf-8") as f: + r = csv.DictReader(f) + for row in r: + ts = datetime.fromisoformat(row["interval_start"].replace("Z", "+00:00")) + px = float(row["sell_raw_price_czk_kwh"]) + out.append((ts, px)) + out.sort(key=lambda x: x[0]) + return out + + +def load_prices_db(date_from: date, date_to: date) -> list[tuple[datetime, float]]: + from datetime import timezone + + try: + import psycopg2 + except ImportError as e: + raise SystemExit("Pro --db instaluj psycopg2-binary nebo použij --price-csv") from e + from zoneinfo import ZoneInfo + + prg = ZoneInfo("Europe/Prague") + t0 = datetime.combine(date_from, datetime.min.time(), tzinfo=prg).astimezone(timezone.utc) + t1 = datetime.combine(date_to, datetime.min.time(), tzinfo=prg).astimezone(timezone.utc) + + conn = psycopg2.connect( + host=os.environ.get("PGHOST", "127.0.0.1"), + port=int(os.environ.get("PGPORT", "5432")), + dbname=os.environ.get("PGDATABASE", "ems"), + user=os.environ.get("PGUSER", os.environ.get("DB_USER", "ems_user")), + password=os.environ.get("PGPASSWORD", os.environ.get("DB_PASSWORD", "")), + ) + cur = conn.cursor() + cur.execute( + """ + SELECT interval_start, sell_raw_price_czk_kwh::float + FROM ems.market_interval_price + WHERE market_source = 'OTE_CZ' + AND interval_start >= %s + AND interval_start < %s + ORDER BY interval_start + """, + (t0, t1), + ) + rows = cur.fetchall() + conn.close() + return [(r[0], float(r[1])) for r in rows] + + +def prices_by_calendar_day( + series: list[tuple[datetime, float]], +) -> dict[date, list[float]]: + """96 hodnot Kč/kWh (raw OTE) na kalendářní den Europe/Prague.""" + from zoneinfo import ZoneInfo + + prg = ZoneInfo("Europe/Prague") + buckets: dict[date, dict[int, float]] = {} + for ts, px in series: + local = ts.astimezone(prg) + d = local.date() + slot = local.hour * 4 + local.minute // 15 + buckets.setdefault(d, {})[slot] = px + out: dict[date, list[float]] = {} + for d, mp in buckets.items(): + if len(mp) < SLOTS_PER_DAY: + continue + out[d] = [mp[i] for i in range(SLOTS_PER_DAY)] + return out + + +def solve_one_day( + pv_wh: Sequence[float], + load_wh: Sequence[float], + p_sell: Sequence[float], + p_buy_flat: float, + e_usable_wh: float, + p_batt_w: float, + site: SiteLimits, + soc_start_wh: float, +) -> tuple[float, float, float, float]: + """ + Vrátí (cash_kc, soc_end_wh, curtailed_wh, discharged_wh_sum). + cash = příjem z exportu − nákup z DS (jen energie, Kč). + """ + e_min = site.soc_min_frac * e_usable_wh + e_max = site.soc_max_frac * e_usable_wh + max_ch = p_batt_w * DT_H + max_dis = p_batt_w * DT_H + max_exp = site.max_export_w * DT_H + max_imp = site.max_import_w * DT_H + + prob = pulp.LpProblem("ems_day", pulp.LpMaximize) + soc = pulp.LpVariable.dicts("soc", range(SLOTS_PER_DAY + 1), lowBound=e_min, upBound=e_max) + ch = pulp.LpVariable.dicts("ch", range(SLOTS_PER_DAY), lowBound=0) + dis = pulp.LpVariable.dicts("dis", range(SLOTS_PER_DAY), lowBound=0) + gexp = pulp.LpVariable.dicts("gexp", range(SLOTS_PER_DAY), lowBound=0) + gimp = pulp.LpVariable.dicts("gimp", range(SLOTS_PER_DAY), lowBound=0) + curt = pulp.LpVariable.dicts("curt", range(SLOTS_PER_DAY), lowBound=0) + + prob += soc[0] == soc_start_wh + + obj = [] + for t in range(SLOTS_PER_DAY): + prob += ch[t] <= max_ch + prob += dis[t] <= max_dis + prob += gexp[t] <= max_exp + prob += gimp[t] <= max_imp + prob += curt[t] <= pv_wh[t] + prob += ( + pv_wh[t] - curt[t] + dis[t] + gimp[t] == load_wh[t] + ch[t] + gexp[t] + ), f"balance_{t}" + prob += ( + soc[t + 1] + == soc[t] + site.eta_charge * ch[t] - dis[t] / site.eta_discharge + ), f"socdyn_{t}" + obj.append(p_sell[t] * gexp[t] / 1000.0 - p_buy_flat * gimp[t] / 1000.0) + + prob += pulp.lpSum(obj) + + solver = pulp.PULP_CBC_CMD(msg=False, timeLimit=60) + prob.solve(solver) + if prob.status != pulp.LpStatusOptimal: + raise RuntimeError(f"LP status {pulp.LpStatus[prob.status]}") + + cash = float(pulp.value(prob.objective)) + soc_end = float(pulp.value(soc[SLOTS_PER_DAY])) + curt_total = sum(float(pulp.value(curt[t])) for t in range(SLOTS_PER_DAY)) + dis_total = sum(float(pulp.value(dis[t])) for t in range(SLOTS_PER_DAY)) + return cash, soc_end, curt_total, dis_total + + +def simulate_year( + days: Iterable[date], + px_day: dict[date, list[float]], + usable_kwh: float, + site: SiteLimits, + sell_margin_fixed: float, + sell_margin_pct: float, + buy_vat_kwh: float, + summer_kwh: float, + winter_kwh: float, + load_kw: float, + shape: Sequence[float], +) -> dict[str, float]: + e_wh = usable_kwh * 1000.0 + p_batt = batt_power_cap_w(usable_kwh, site) + load_wh = daily_load_wh(load_kw) + cash_total = 0.0 + curt_total = 0.0 + dis_total = 0.0 + soc_state = 0.5 * (site.soc_min_frac + site.soc_max_frac) * e_wh + n_days = 0 + for d in days: + if d not in px_day: + continue + raw = px_day[d] + p_sell = [effective_sell_kc_kwh(x, sell_margin_fixed, sell_margin_pct) for x in raw] + pv_wh = daily_pv_wh(d, summer_kwh, winter_kwh, shape) + cash, soc_state, curt, dis = solve_one_day( + pv_wh, load_wh, p_sell, buy_vat_kwh, e_wh, p_batt, site, soc_state + ) + cash_total += cash + curt_total += curt + dis_total += dis + n_days += 1 + feq = (dis_total / e_wh / n_days) if n_days and e_wh > 0 else 0.0 + return { + "cash_kc": cash_total, + "days": float(n_days), + "curt_wh": curt_total, + "dis_wh": dis_total, + "feq_cycles_per_day": feq, + } + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--db", action="store_true", help="Načti OTE z Postgres (env PG* / DB_*)") + ap.add_argument("--price-csv", type=str, default="", help="CSV: interval_start, sell_raw_price_czk_kwh") + ap.add_argument("--date-from", type=str, required=True) + ap.add_argument("--date-to", type=str, required=True) + ap.add_argument("--battery-kwh", type=float, nargs="+", required=True, help="Užitkové kWh (např. 12.5 32 48)") + ap.add_argument("--load-kw", type=float, default=1.0, help="Průměrný odběr (konstanta přes den)") + ap.add_argument("--pv-daily-kwh-summer", type=float, default=50.0) + ap.add_argument("--pv-daily-kwh-winter", type=float, default=10.0) + ap.add_argument("--sell-margin-fixed", type=float, default=-0.02) + ap.add_argument("--sell-margin-pct", type=float, default=0.0) + ap.add_argument("--buy-vat-kwh", type=float, default=4.443, help="Efektivní nákup Kč/kWh vč. DPH (flat screening)") + ap.add_argument("--max-export-w", type=float, default=16_000.0) + ap.add_argument("--max-import-w", type=float, default=17_000.0) + ap.add_argument("--inv-batt-max-w", type=float, default=12_000.0) + ap.add_argument("--c-rate", type=float, default=0.5) + ap.add_argument("--capex-per-kwh", type=float, default=0.0, help="CAPEX za 1 kWh rozšíření; vypíše jednoduchou návratnost vs. nejmenší baterie") + args = ap.parse_args() + + d0 = date.fromisoformat(args.date_from) + d1 = date.fromisoformat(args.date_to) + if args.db: + series = load_prices_db(d0, d1) + elif args.price_csv: + series = load_prices_csv(args.price_csv) + else: + ap.error("Zadej --db nebo --price-csv") + + px_day = prices_by_calendar_day(series) + shape = pv_shape_96() + site = SiteLimits( + max_export_w=args.max_export_w, + max_import_w=args.max_import_w, + inv_batt_max_w=args.inv_batt_max_w, + c_rate=args.c_rate, + ) + + day_list = [d0 + timedelta(days=i) for i in range((d1 - d0).days)] + + results = [] + for kwh in sorted(args.battery_kwh): + r = simulate_year( + day_list, + px_day, + kwh, + site, + args.sell_margin_fixed, + args.sell_margin_pct, + args.buy_vat_kwh, + args.pv_daily_kwh_summer, + args.pv_daily_kwh_winter, + args.load_kw, + shape, + ) + results.append((kwh, r)) + + baseline_kwh = min(args.battery_kwh) + base = dict(results)[baseline_kwh] + + print("Parametry: prodej = OTE + sell_margin_fixed (+ %), nákup = flat buy_vat_kwh") + print(f" FVE tvar = syntetický den, léto {args.pv_daily_kwh_summer} kWh/d, zima {args.pv_daily_kwh_winter} kWh/d, load {args.load_kw} kW") + print(f" Limity: export {args.max_export_w} W, import {args.max_import_w} W, P_batt = min({args.c_rate}*E_kWh, {args.inv_batt_max_w} W)") + print() + + print(f"{'kWh':>8} {'P_batt_kW':>10} {'cash_kc/rok':>14} {'Δ vs min':>12} {'curt_MWh/y':>12} {'Feq/den':>8}") + for kwh, r in results: + pkw = batt_power_cap_w(kwh, site) / 1000.0 + days = max(int(r["days"]), 1) + cash_y = r["cash_kc"] * (365.0 / days) + curt_mwh = r["curt_wh"] / 1e6 * (365.0 / days) + delta = cash_y - base["cash_kc"] * (365.0 / days) if kwh != baseline_kwh else 0.0 + print( + f"{kwh:8.1f} {pkw:10.2f} {cash_y:14.0f} {delta:12.0f} {curt_mwh:12.2f} {r['feq_cycles_per_day']:8.2f}" + ) + + if args.capex_per_kwh > 0: + print() + base_cash = base["cash_kc"] * (365.0 / max(int(base["days"]), 1)) + for kwh, r in results: + if kwh <= baseline_kwh: + continue + cash_y = r["cash_kc"] * (365.0 / max(int(r["days"]), 1)) + delta = cash_y - base_cash + extra_kwh = kwh - baseline_kwh + capex = extra_kwh * args.capex_per_kwh + if delta > 0: + years = capex / delta + print( + f"vs {baseline_kwh} kWh → +{extra_kwh:.0f} kWh CAPEX ~{capex:,.0f} Kč, " + f"odhad +{delta:,.0f} Kč/rok → návratnost ~{years:.1f} r" + ) + else: + print(f"vs {baseline_kwh} kWh → +{extra_kwh:.0f} kWh: model neukazuje vyšší roční cash ({delta:,.0f} Kč/rok)") + + +if __name__ == "__main__": + main() diff --git a/scripts/analysis/data/ote_2026-04-09_11_12.csv b/scripts/analysis/data/ote_2026-04-09_11_12.csv new file mode 100644 index 0000000..f2cf5ef --- /dev/null +++ b/scripts/analysis/data/ote_2026-04-09_11_12.csv @@ -0,0 +1,288 @@ +2026-04-09,00:00,2.6720 +2026-04-09,00:15,2.5505 +2026-04-09,00:30,2.4905 +2026-04-09,00:45,2.4618 +2026-04-09,01:00,2.4598 +2026-04-09,01:15,2.4830 +2026-04-09,01:30,2.4660 +2026-04-09,01:45,2.4603 +2026-04-09,02:00,2.5270 +2026-04-09,02:15,2.4995 +2026-04-09,02:30,2.5213 +2026-04-09,02:45,2.5305 +2026-04-09,03:00,2.5265 +2026-04-09,03:15,2.5483 +2026-04-09,03:30,2.5655 +2026-04-09,03:45,2.6203 +2026-04-09,04:00,2.5300 +2026-04-09,04:15,2.6015 +2026-04-09,04:30,2.6713 +2026-04-09,04:45,2.8440 +2026-04-09,05:00,2.6540 +2026-04-09,05:15,2.7750 +2026-04-09,05:30,3.0300 +2026-04-09,05:45,3.2593 +2026-04-09,06:00,2.9765 +2026-04-09,06:15,3.4128 +2026-04-09,06:30,3.7728 +2026-04-09,06:45,3.8188 +2026-04-09,07:00,3.9350 +2026-04-09,07:15,3.8950 +2026-04-09,07:30,3.6658 +2026-04-09,07:45,3.2233 +2026-04-09,08:00,3.9688 +2026-04-09,08:15,3.3950 +2026-04-09,08:30,2.9795 +2026-04-09,08:45,2.4095 +2026-04-09,09:00,3.1465 +2026-04-09,09:15,2.6005 +2026-04-09,09:30,2.1475 +2026-04-09,09:45,1.6428 +2026-04-09,10:00,2.0153 +2026-04-09,10:15,1.6153 +2026-04-09,10:30,0.9560 +2026-04-09,10:45,0.2468 +2026-04-09,11:00,0.1235 +2026-04-09,11:15,0.0003 +2026-04-09,11:30,0.0000 +2026-04-09,11:45,-0.0003 +2026-04-09,12:00,-0.0003 +2026-04-09,12:15,-0.0025 +2026-04-09,12:30,-0.0313 +2026-04-09,12:45,-0.0760 +2026-04-09,13:00,-0.0808 +2026-04-09,13:15,-0.1385 +2026-04-09,13:30,-0.2323 +2026-04-09,13:45,-0.3113 +2026-04-09,14:00,-0.2458 +2026-04-09,14:15,-0.2328 +2026-04-09,14:30,-0.1720 +2026-04-09,14:45,-0.1265 +2026-04-09,15:00,-0.0748 +2026-04-09,15:15,-0.0053 +2026-04-09,15:30,0.0035 +2026-04-09,15:45,0.0003 +2026-04-09,16:00,0.0003 +2026-04-09,16:15,0.1078 +2026-04-09,16:30,0.9075 +2026-04-09,16:45,1.7903 +2026-04-09,17:00,0.7188 +2026-04-09,17:15,1.9618 +2026-04-09,17:30,2.6005 +2026-04-09,17:45,3.1290 +2026-04-09,18:00,2.4015 +2026-04-09,18:15,2.8648 +2026-04-09,18:30,3.1983 +2026-04-09,18:45,3.4253 +2026-04-09,19:00,3.4435 +2026-04-09,19:15,3.6148 +2026-04-09,19:30,3.8523 +2026-04-09,19:45,3.8893 +2026-04-09,20:00,4.1890 +2026-04-09,20:15,3.8510 +2026-04-09,20:30,3.5613 +2026-04-09,20:45,3.1660 +2026-04-09,21:00,3.3128 +2026-04-09,21:15,3.0373 +2026-04-09,21:30,2.6313 +2026-04-09,21:45,2.2740 +2026-04-09,22:00,2.9655 +2026-04-09,22:15,2.7910 +2026-04-09,22:30,2.5323 +2026-04-09,22:45,2.2793 +2026-04-09,23:00,2.4410 +2026-04-09,23:15,2.2813 +2026-04-09,23:30,2.1535 +2026-04-09,23:45,1.9268 +2026-04-11,00:00,3.5698 +2026-04-11,00:15,3.2533 +2026-04-11,00:30,2.9323 +2026-04-11,00:45,2.7985 +2026-04-11,01:00,3.2713 +2026-04-11,01:15,3.1870 +2026-04-11,01:30,2.9530 +2026-04-11,01:45,2.9568 +2026-04-11,02:00,2.9888 +2026-04-11,02:15,2.9133 +2026-04-11,02:30,2.9605 +2026-04-11,02:45,2.9473 +2026-04-11,03:00,3.0010 +2026-04-11,03:15,2.9330 +2026-04-11,03:30,2.9580 +2026-04-11,03:45,2.8618 +2026-04-11,04:00,2.9760 +2026-04-11,04:15,2.9570 +2026-04-11,04:30,2.9235 +2026-04-11,04:45,2.8015 +2026-04-11,05:00,2.9408 +2026-04-11,05:15,2.8243 +2026-04-11,05:30,2.8233 +2026-04-11,05:45,2.7850 +2026-04-11,06:00,2.8100 +2026-04-11,06:15,2.8503 +2026-04-11,06:30,2.8383 +2026-04-11,06:45,2.6828 +2026-04-11,07:00,2.9005 +2026-04-11,07:15,2.8145 +2026-04-11,07:30,2.4048 +2026-04-11,07:45,1.7628 +2026-04-11,08:00,2.7298 +2026-04-11,08:15,2.3223 +2026-04-11,08:30,1.6730 +2026-04-11,08:45,0.5990 +2026-04-11,09:00,1.4445 +2026-04-11,09:15,0.5615 +2026-04-11,09:30,0.2000 +2026-04-11,09:45,0.0003 +2026-04-11,10:00,0.1933 +2026-04-11,10:15,0.0000 +2026-04-11,10:30,-0.0305 +2026-04-11,10:45,-0.2223 +2026-04-11,11:00,-0.0308 +2026-04-11,11:15,0.0000 +2026-04-11,11:30,0.0250 +2026-04-11,11:45,0.1450 +2026-04-11,12:00,-0.4685 +2026-04-11,12:15,-1.0418 +2026-04-11,12:30,-0.9168 +2026-04-11,12:45,-1.2465 +2026-04-11,13:00,-1.5560 +2026-04-11,13:15,-1.7465 +2026-04-11,13:30,-2.1328 +2026-04-11,13:45,-2.4545 +2026-04-11,14:00,-1.8465 +2026-04-11,14:15,-1.8458 +2026-04-11,14:30,-1.7578 +2026-04-11,14:45,-1.4893 +2026-04-11,15:00,-1.7045 +2026-04-11,15:15,-1.3535 +2026-04-11,15:30,-1.0890 +2026-04-11,15:45,-0.6923 +2026-04-11,16:00,-0.7728 +2026-04-11,16:15,-0.3718 +2026-04-11,16:30,-0.2075 +2026-04-11,16:45,-0.3675 +2026-04-11,17:00,0.3510 +2026-04-11,17:15,-0.1903 +2026-04-11,17:30,-0.0250 +2026-04-11,17:45,0.1198 +2026-04-11,18:00,-0.4573 +2026-04-11,18:15,1.8173 +2026-04-11,18:30,2.3750 +2026-04-11,18:45,2.9240 +2026-04-11,19:00,2.3338 +2026-04-11,19:15,2.6563 +2026-04-11,19:30,2.8090 +2026-04-11,19:45,3.0140 +2026-04-11,20:00,2.8350 +2026-04-11,20:15,2.4450 +2026-04-11,20:30,2.4565 +2026-04-11,20:45,2.2635 +2026-04-11,21:00,3.0835 +2026-04-11,21:15,2.1543 +2026-04-11,21:30,1.9325 +2026-04-11,21:45,1.8295 +2026-04-11,22:00,3.3630 +2026-04-11,22:15,2.5945 +2026-04-11,22:30,1.7628 +2026-04-11,22:45,1.4008 +2026-04-11,23:00,2.4143 +2026-04-11,23:15,1.7565 +2026-04-11,23:30,1.5728 +2026-04-11,23:45,1.5035 +2026-04-12,00:00,1.5390 +2026-04-12,00:15,1.5708 +2026-04-12,00:30,1.7003 +2026-04-12,00:45,1.7023 +2026-04-12,01:00,1.7770 +2026-04-12,01:15,1.7503 +2026-04-12,01:30,1.8833 +2026-04-12,01:45,2.0250 +2026-04-12,02:00,1.8593 +2026-04-12,02:15,1.9398 +2026-04-12,02:30,1.9483 +2026-04-12,02:45,2.0395 +2026-04-12,03:00,2.0418 +2026-04-12,03:15,2.1465 +2026-04-12,03:30,2.1558 +2026-04-12,03:45,2.2473 +2026-04-12,04:00,2.3520 +2026-04-12,04:15,2.4003 +2026-04-12,04:30,2.4923 +2026-04-12,04:45,2.5645 +2026-04-12,05:00,2.4895 +2026-04-12,05:15,2.5743 +2026-04-12,05:30,2.5745 +2026-04-12,05:45,2.5400 +2026-04-12,06:00,2.6265 +2026-04-12,06:15,2.6030 +2026-04-12,06:30,2.5383 +2026-04-12,06:45,2.4885 +2026-04-12,07:00,2.7490 +2026-04-12,07:15,2.5110 +2026-04-12,07:30,2.3325 +2026-04-12,07:45,2.2773 +2026-04-12,08:00,2.1050 +2026-04-12,08:15,1.8573 +2026-04-12,08:30,1.7498 +2026-04-12,08:45,1.3183 +2026-04-12,09:00,1.4575 +2026-04-12,09:15,1.3108 +2026-04-12,09:30,0.8415 +2026-04-12,09:45,0.6773 +2026-04-12,10:00,0.8375 +2026-04-12,10:15,0.4075 +2026-04-12,10:30,0.3775 +2026-04-12,10:45,0.3193 +2026-04-12,11:00,0.4165 +2026-04-12,11:15,0.3638 +2026-04-12,11:30,0.3085 +2026-04-12,11:45,0.2470 +2026-04-12,12:00,0.2023 +2026-04-12,12:15,0.1435 +2026-04-12,12:30,0.1298 +2026-04-12,12:45,0.0615 +2026-04-12,13:00,0.0473 +2026-04-12,13:15,0.0173 +2026-04-12,13:30,0.0015 +2026-04-12,13:45,0.0013 +2026-04-12,14:00,0.0000 +2026-04-12,14:15,0.0278 +2026-04-12,14:30,0.0688 +2026-04-12,14:45,0.0628 +2026-04-12,15:00,0.1483 +2026-04-12,15:15,0.2505 +2026-04-12,15:30,0.4290 +2026-04-12,15:45,0.7538 +2026-04-12,16:00,0.8148 +2026-04-12,16:15,1.4820 +2026-04-12,16:30,1.5620 +2026-04-12,16:45,2.4738 +2026-04-12,17:00,1.5428 +2026-04-12,17:15,2.2933 +2026-04-12,17:30,3.1040 +2026-04-12,17:45,3.4375 +2026-04-12,18:00,2.8968 +2026-04-12,18:15,3.1310 +2026-04-12,18:30,3.4105 +2026-04-12,18:45,3.6880 +2026-04-12,19:00,3.3903 +2026-04-12,19:15,3.3798 +2026-04-12,19:30,3.6125 +2026-04-12,19:45,3.9263 +2026-04-12,20:00,3.7030 +2026-04-12,20:15,3.6903 +2026-04-12,20:30,3.6473 +2026-04-12,20:45,3.5030 +2026-04-12,21:00,3.5230 +2026-04-12,21:15,3.3825 +2026-04-12,21:30,3.4325 +2026-04-12,21:45,3.3353 +2026-04-12,22:00,3.4898 +2026-04-12,22:15,3.3675 +2026-04-12,22:30,3.2180 +2026-04-12,22:45,2.9855 +2026-04-12,23:00,3.3763 +2026-04-12,23:15,2.9848 +2026-04-12,23:30,2.9855 +2026-04-12,23:45,2.7558 diff --git a/scripts/analysis/join_inverter_ote_snapshot.py b/scripts/analysis/join_inverter_ote_snapshot.py new file mode 100644 index 0000000..910f303 --- /dev/null +++ b/scripts/analysis/join_inverter_ote_snapshot.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +"""Join Deye inverter export (wide xlsx) with OTE 15min sell prices for BA81-style analysis. + +OTE CSV: regenerate from EMS DB (MCP or psql), example: + + SELECT string_agg( + to_char((interval_start AT TIME ZONE 'Europe/Prague')::date, 'YYYY-MM-DD') || ',' || + to_char(interval_start AT TIME ZONE 'Europe/Prague', 'HH24:MI') || ',' || + trim(to_char(sell_raw_price_czk_kwh, 'FM9999990.0000')), + chr(10) ORDER BY interval_start + ) + FROM ems.market_interval_price + WHERE market_source = 'OTE_CZ' + AND (interval_start AT TIME ZONE 'Europe/Prague')::date IN (...); + +Convention in sample logs: negative Battery Power(W) ≈ charging, positive ≈ discharging. +Total Grid Power(W): small positive ≈ little/no export (sign per site firmware). + +Requires: openpyxl. Use read_only=False (these exports report max_row=1 in read_only mode). +""" +from __future__ import annotations + +import argparse +import statistics as st +from collections import defaultdict +from datetime import datetime +from pathlib import Path + +import openpyxl + +COLS = [ + "Time", + "Total Solar Power(W)", + "Total Inverter Output Power(W)", + "Total Grid Power(W)", + "Battery Power(W)", + "SoC(%)", +] + + +def load_ote_csv(path: Path) -> dict[tuple[str, str], float]: + ote: dict[tuple[str, str], float] = {} + for line in path.read_text().splitlines(): + line = line.strip() + if not line: + continue + d, hm, s = line.split(",") + ote[(d, hm)] = float(s) + return ote + + +def floor_15(dt: datetime) -> datetime: + m = (dt.minute // 15) * 15 + return dt.replace(minute=m, second=0, microsecond=0) + + +def slot_key(dt: datetime) -> tuple[str, str]: + f = floor_15(dt) + return f.strftime("%Y-%m-%d"), f.strftime("%H:%M") + + +def load_inverter_rows(fp: Path) -> list[dict[str, object]]: + wb = openpyxl.load_workbook(fp, read_only=False, data_only=True) + ws = wb.active + it = ws.iter_rows(values_only=True) + header = next(it) + idx = {str(h).strip(): i for i, h in enumerate(header) if h} + rows: list[dict[str, object]] = [] + for r in it: + if not r or r[idx["Time"]] is None: + continue + rows.append({c: r[idx[c]] for c in COLS}) + wb.close() + return rows + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--ote-csv", type=Path, required=True) + p.add_argument("xlsx", type=Path, nargs="+") + args = p.parse_args() + ote = load_ote_csv(args.ote_csv) + + for fp in args.xlsx: + data = load_inverter_rows(fp) + neg: list[tuple[float, datetime, dict]] = [] + for r in data: + t = r["Time"] + if isinstance(t, str): + t = datetime.strptime(t, "%Y/%m/%d %H:%M:%S") + dk, hm = slot_key(t) + sell = ote.get((dk, hm)) + if sell is None or sell >= 0: + continue + neg.append((sell, t, r)) + + print(f"\n=== {fp.name} rows={len(data)} OTE sell<0 samples={len(neg)}") + if not neg: + continue + socs = [float(x[2]["SoC(%)"]) for x in neg] + grids = [float(x[2]["Total Grid Power(W)"]) for x in neg] + bats = [float(x[2]["Battery Power(W)"]) for x in neg] + sols = [float(x[2]["Total Solar Power(W)"]) for x in neg] + print(f" SoC %: mean={st.mean(socs):.1f} min={min(socs):.0f} max={max(socs):.0f}") + print(f" Grid W: mean={st.mean(grids):.0f} med={st.median(grids):.0f}") + print(f" Bat W: mean={st.mean(bats):.0f} med={st.median(bats):.0f}") + print(f" Solar W: mean={st.mean(sols):.0f} med={st.median(sols):.0f}") + + buckets: dict[str, list] = defaultdict(list) + for sell, t, r in neg: + if t.hour < 9 or t.hour > 18: + continue + _, hm = slot_key(t) + buckets[hm].append((sell, r)) + print(" 15min buckets (OTE<0, 09-18h) medians:") + for hm in sorted(buckets.keys()): + b = buckets[hm] + sell = b[0][0] + socs_b = [float(x[1]["SoC(%)"]) for x in b] + print( + f" {hm} sell={sell:+.3f} n={len(b):2d} " + f"SoC_med={st.median(socs_b):.0f}% " + f"Pgrid_med={st.median([float(x[1]['Total Grid Power(W)']) for x in b]):.0f}W " + f"Psol_med={st.median([float(x[1]['Total Solar Power(W)']) for x in b]):.0f}W" + ) + + +if __name__ == "__main__": + main()