tune microcycling

scripts/analysis/battery_sizing_screen.py

@@ -13,22 +13,27 @@ Připojení k DB (deploy / Docker):
   - Nebo ``DATABASE_URL`` / ``postgresql://USER:PASS@HOST:5432/ems`` (na hostu HOST=127.0.0.1
     nebo EMS_DB_BIND, ne ``db`` — to je jen uvnitř Docker sítě).
 
-Příklad:
-  python3 scripts/analysis/battery_sizing_screen.py \\
-    --db \\
+Příklad (syntetická FVE, flat nákup):
+  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 \\
+    --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
+    --sell-margin-fixed -0.02 --buy-vat-kwh 4.443 --capex-per-kwh 9000
+
+Příklad (PVGIS měsíční E_d + NT/VT):
+  python3 scripts/analysis/battery_sizing_screen.py --db \\
+    --pvgis-csv pole_A.csv --pvgis-csv pole_B.csv \\
+    --buy-nt-kwh 5.25 --buy-vt-surcharge-kwh 2.0 --nt-from-hour 22 --nt-to-hour 6 \\
+    ... (ostatní jako výše)
 
 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.
+Omezení modelu: FVE buď syntetický denní tvar (--pv-daily-kwh-*), nebo součet měsíčních
+E_d z PVGIS CSV (--pvgis-csv, opakovat pro více orientací); denní energie = E_d měsíce
+× normalizovaný tvar (stejný profil každý den v měsíci). Nákup: buď flat (--buy-vat-kwh),
+nebo NT/VT podle hodin Europe/Prague: --buy-nt-kwh, VT = NT + --buy-vt-surcharge-kwh,
+okno NT --nt-from-hour až --nt-to-hour (přes půlnoc, pokud from > to). Mikroinvertory / GEN
+nejsou; zelený bonus není v účelové funkci. Výsledek = screening, ne nabídka.
 """
 from __future__ import annotations
 
@@ -40,7 +45,7 @@ import sys
 from dataclasses import dataclass
 from datetime import date, datetime, timedelta
 from pathlib import Path
-from typing import Iterable, Sequence
+from typing import Iterable, Sequence, Mapping
 
 try:
     import pulp
@@ -93,6 +98,72 @@ def daily_pv_wh(d: date, summer_kwh: float, winter_kwh: float, shape: Sequence[f
     return [base * 1000.0 * sh for sh in shape]
 
 
+def load_pvgis_monthly_ed_kwh(path: Path) -> dict[int, float]:
+    """Z PVGIS CSV (Fixed angle) načte E_d [kWh/d] pro měsíce 1–12."""
+    text = path.read_text(encoding="utf-8", errors="replace").splitlines()
+    start: int | None = None
+    for i, line in enumerate(text):
+        if line.strip().startswith("Fixed angle"):
+            start = i + 2
+            break
+    if start is None:
+        raise ValueError(f"PVGIS: řádek 'Fixed angle' nenalezen: {path}")
+    out: dict[int, float] = {}
+    for line in text[start:]:
+        cells = [c.strip() for c in line.split("\t") if c.strip() != ""]
+        if not cells:
+            continue
+        if cells[0] == "Year":
+            break
+        try:
+            month = int(cells[0])
+        except ValueError:
+            continue
+        if not (1 <= month <= 12):
+            continue
+        out[month] = float(cells[1].replace(",", "."))
+    if len(out) != 12:
+        raise ValueError(f"PVGIS: očekáváno 12 měsíců E_d v {path}, mám {sorted(out.keys())}")
+    return out
+
+
+def merge_pvgis_monthly_ed_kwh(paths: Sequence[Path]) -> dict[int, float]:
+    """Sečte E_d jednotlivých polí (např. dvě orientace)."""
+    total = {m: 0.0 for m in range(1, 13)}
+    for p in paths:
+        part = load_pvgis_monthly_ed_kwh(Path(p))
+        for m in range(1, 13):
+            total[m] += part[m]
+    return total
+
+
+def daily_pv_wh_monthly(d: date, monthly_ed_kwh: Mapping[int, float], shape: Sequence[float]) -> list[float]:
+    kwh = float(monthly_ed_kwh[d.month])
+    return [kwh * 1000.0 * sh for sh in shape]
+
+
+def buy_prices_96_nt_vt(
+    nt_kwh: float,
+    vt_kwh: float,
+    nt_from_hour: int,
+    nt_to_hour: int,
+) -> list[float]:
+    """
+    96 cen nákupu [Kč/kWh] podle začátku 15min slotu (hodina 0–23, Europe/Prague).
+    Pokud nt_from_hour > nt_to_hour: NT pro hodiny >= from nebo < to (přes půlnoc).
+    Jinak NT pro from <= h < to.
+    """
+    out: list[float] = []
+    for t in range(SLOTS_PER_DAY):
+        h = t // 4
+        if nt_from_hour > nt_to_hour:
+            is_nt = h >= nt_from_hour or h < nt_to_hour
+        else:
+            is_nt = nt_from_hour <= h < nt_to_hour
+        out.append(nt_kwh if is_nt else vt_kwh)
+    return out
+
+
 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
@@ -221,7 +292,7 @@ def solve_one_day(
     pv_wh: Sequence[float],
     load_wh: Sequence[float],
     p_sell: Sequence[float],
-    p_buy_flat: float,
+    p_buy: Sequence[float],
     e_usable_wh: float,
     p_batt_w: float,
     site: SiteLimits,
@@ -262,7 +333,7 @@ def solve_one_day(
             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)
+        obj.append(p_sell[t] * gexp[t] / 1000.0 - p_buy[t] * gimp[t] / 1000.0)
 
     prob += pulp.lpSum(obj)
 
@@ -285,15 +356,23 @@ def simulate_year(
     site: SiteLimits,
     sell_margin_fixed: float,
     sell_margin_pct: float,
-    buy_vat_kwh: float,
+    buy_flat_kwh: float,
+    buy_prices_96: Sequence[float] | None,
     summer_kwh: float,
     winter_kwh: float,
     load_kw: float,
     shape: Sequence[float],
+    monthly_ed_kwh: Mapping[int, float] | None,
 ) -> 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)
+    if buy_prices_96 is not None:
+        if len(buy_prices_96) != SLOTS_PER_DAY:
+            raise ValueError("buy_prices_96 musí mít 96 hodnot")
+        p_buy_day: Sequence[float] = buy_prices_96
+    else:
+        p_buy_day = [buy_flat_kwh] * SLOTS_PER_DAY
     cash_total = 0.0
     curt_total = 0.0
     dis_total = 0.0
@@ -304,9 +383,12 @@ def simulate_year(
             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)
+        if monthly_ed_kwh is not None:
+            pv_wh = daily_pv_wh_monthly(d, monthly_ed_kwh, shape)
+        else:
+            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
+            pv_wh, load_wh, p_sell, p_buy_day, e_wh, p_batt, site, soc_state
         )
         cash_total += cash
         curt_total += curt
@@ -346,7 +428,33 @@ def main() -> None:
     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(
+        "--buy-vat-kwh",
+        type=float,
+        default=4.443,
+        help="Flat nákup Kč/kWh (když není --buy-nt-kwh)",
+    )
+    ap.add_argument(
+        "--buy-nt-kwh",
+        type=float,
+        default=None,
+        help="NT cena Kč/kWh; VT = NT + --buy-vt-surcharge-kwh; okno --nt-from-hour / --nt-to-hour (Europe/Prague)",
+    )
+    ap.add_argument(
+        "--buy-vt-surcharge-kwh",
+        type=float,
+        default=0.0,
+        help="Příplatek VT oproti NT (jako buy_fixed_vt_surcharge v EMS)",
+    )
+    ap.add_argument("--nt-from-hour", type=int, default=22, help="Začátek NT (hodina 0–23)")
+    ap.add_argument("--nt-to-hour", type=int, default=6, help="Konec NT: první hodina VT (0–23); přes půlnoc pokud from > to")
+    ap.add_argument(
+        "--pvgis-csv",
+        action="append",
+        default=[],
+        metavar="PATH",
+        help="PVGIS měsíční E_d (Fixed angle); opakovat pro více polí/orientací, energie se sečte",
+    )
     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)
@@ -379,6 +487,23 @@ def main() -> None:
         c_rate=args.c_rate,
     )
 
+    monthly_ed: dict[int, float] | None = None
+    if args.pvgis_csv:
+        monthly_ed = merge_pvgis_monthly_ed_kwh([Path(p) for p in args.pvgis_csv])
+
+    if args.buy_nt_kwh is not None:
+        vt = args.buy_nt_kwh + args.buy_vt_surcharge_kwh
+        buy_prices_96 = buy_prices_96_nt_vt(
+            args.buy_nt_kwh,
+            vt,
+            args.nt_from_hour,
+            args.nt_to_hour,
+        )
+        buy_flat = args.buy_vat_kwh
+    else:
+        buy_prices_96 = None
+        buy_flat = args.buy_vat_kwh
+
     day_list = [d0 + timedelta(days=i) for i in range((d1 - d0).days)]
 
     results = []
@@ -390,19 +515,40 @@ def main() -> None:
             site,
             args.sell_margin_fixed,
             args.sell_margin_pct,
-            args.buy_vat_kwh,
+            buy_flat,
+            buy_prices_96,
             args.pv_daily_kwh_summer,
             args.pv_daily_kwh_winter,
             args.load_kw,
             shape,
+            monthly_ed,
         )
         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("Parametry: prodej = OTE + sell_margin_fixed (+ %)")
+    if buy_prices_96 is not None:
+        vt_show = args.buy_nt_kwh + args.buy_vt_surcharge_kwh
+        print(
+            f"  Nákup = NT/VT: NT {args.buy_nt_kwh} Kč/kWh, VT {vt_show} Kč/kWh "
+            f"(okno NT {args.nt_from_hour:02d}–{args.nt_to_hour:02d} h lokální)"
+        )
+    else:
+        print(f"  Nákup = flat {args.buy_vat_kwh} Kč/kWh")
+    if monthly_ed is not None:
+        edv = [monthly_ed[m] for m in range(1, 13)]
+        print(
+            f"  FVE = PVGIS měsíční E_d (součet {len(args.pvgis_csv)} souborů), "
+            f"rozsah {min(edv):.1f}–{max(edv):.1f} kWh/d, denní tvar = syntetika"
+        )
+    else:
+        print(
+            f"  FVE = syntetický den, léto {args.pv_daily_kwh_summer} kWh/d, "
+            f"zima {args.pv_daily_kwh_winter} kWh/d"
+        )
+    print(f"  Load (konstanta) {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()