904 lines
33 KiB
Python
904 lines
33 KiB
Python
#!/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řipojení k DB (deploy / Docker):
|
||
- Postgres v compose poslouchá na ``EMS_DB_BIND:5432`` (výchozí 127.0.0.1). ``connection refused``
|
||
= služba ``db`` neběží, nebo je port vázaný jen na jinou IP (WireGuard) → nastav stejný host.
|
||
- Skript načte první nalezené ``.env`` z ``…/ems-deploy/.env`` nebo ``…/app/.env`` (není-li
|
||
``--no-auto-env``) a doplní ``PGUSER``/``PGPASSWORD`` z ``DB_USER``/``DB_PASSWORD``.
|
||
- 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 (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 \\
|
||
--pv-daily-kwh-summer 55 --pv-daily-kwh-winter 12 \\
|
||
--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)
|
||
|
||
Příklad (čistá arbitráž, nákup = spot + fixní adder + distribuce/poplatky):
|
||
python3 scripts/analysis/battery_sizing_screen.py --db \\
|
||
--load-kw 0 --pv-daily-kwh-summer 0 --pv-daily-kwh-winter 0 \\
|
||
--buy-spot-add-fixed-kwh 0.25 \\
|
||
--buy-distribution-kwh 1.80 --buy-other-fees-kwh 0.20 --buy-vat-multiplier 1.21 \\
|
||
... (ostatní jako výše)
|
||
|
||
Příklad (nákup jako home-01, prodej spot −0,30 dle živé DB):
|
||
python3 scripts/analysis/battery_sizing_screen.py --db \\
|
||
--buy-home-01 --sell-margin-fixed -0.30 ... (ostatní jako výše)
|
||
|
||
Vyžaduje: pip install pulp (volitelně psycopg2 pro --db).
|
||
|
||
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: flat (--buy-vat-kwh),
|
||
NT/VT podle hodin Europe/Prague (--buy-nt-kwh, VT = NT + --buy-vt-surcharge-kwh),
|
||
nebo od raw OTE spotu: --buy-spot-add-fixed-kwh / --buy-spot-asym-pct; u všech režimů
|
||
lze přičíst --buy-distribution-kwh a --buy-other-fees-kwh a výslednou cenu násobit
|
||
--buy-vat-multiplier. Model explicitně zakazuje současný import+export a současné
|
||
nabíjení+vybíjení v jednom slotu. Dlouhé běhy MILP lze řídit přes
|
||
--solver-time-limit-sec a průběžný tisk přes --progress-every-days. Mikroinvertory /
|
||
GEN nejsou; zelený bonus není v účelové funkci. 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 pathlib import Path
|
||
from time import perf_counter
|
||
from typing import Iterable, Sequence, Mapping
|
||
|
||
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
|
||
|
||
|
||
# home-01 (živá DB site_market_config + V016 tarif/HDO): ověř MCP před změnou
|
||
HOME01_BUY_MARGIN_FIXED_KWH = 0.0
|
||
HOME01_BUY_MARGIN_PCT = 9.0
|
||
HOME01_SYSTEM_SERVICES_KWH = 0.192
|
||
HOME01_OTE_FEE_KWH = 0.001
|
||
HOME01_DIST_NT_KWH = 0.2243
|
||
HOME01_DIST_VT_KWH = 0.74987
|
||
HOME01_VAT_MULTIPLIER = 1.21
|
||
HOME01_SELL_MARGIN_FIXED_KWH = -0.30
|
||
# VT okna Europe/Prague (každý den): 09–10, 12–13, 16–17, 20–21
|
||
HOME01_VT_SLOT_MINUTES: tuple[tuple[int, int], ...] = (
|
||
(9 * 60, 10 * 60),
|
||
(12 * 60, 13 * 60),
|
||
(16 * 60, 17 * 60),
|
||
(20 * 60, 21 * 60),
|
||
)
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class BuyPricingConfig:
|
||
mode: str = "flat"
|
||
flat_kwh: float = 4.443
|
||
nt_kwh: float | None = None
|
||
vt_kwh: float | None = None
|
||
nt_from_hour: int = 22
|
||
nt_to_hour: int = 6
|
||
spot_add_fixed_kwh: float | None = None
|
||
spot_asym_pct: float | None = None
|
||
distribution_kwh: float = 0.0
|
||
other_fees_kwh: float = 0.0
|
||
vat_multiplier: float = 1.0
|
||
|
||
|
||
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 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
|
||
|
||
|
||
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 effective_buy_spot_add_fixed_kc_kwh(raw_ote: float, add_fixed_kwh: float) -> float:
|
||
return raw_ote + add_fixed_kwh
|
||
|
||
|
||
def effective_buy_spot_asym_pct_kc_kwh(raw_ote: float, asym_pct: float) -> float:
|
||
if raw_ote >= 0:
|
||
return raw_ote * (1.0 + asym_pct / 100.0)
|
||
return raw_ote * (1.0 - asym_pct / 100.0)
|
||
|
||
|
||
def slot_is_vt_home01(slot_index: int) -> bool:
|
||
"""15min slot 0..95 od půlnoci (Europe/Prague) — VT okna jako u home-01 HDO."""
|
||
mins = (slot_index // 4) * 60 + (slot_index % 4) * 15
|
||
for start_m, end_m in HOME01_VT_SLOT_MINUTES:
|
||
if start_m <= mins < end_m:
|
||
return True
|
||
return False
|
||
|
||
|
||
def effective_buy_home01_kc_kwh(raw_ote: float, is_vt: bool) -> float:
|
||
"""Stejná struktura jako ems.fn_effective_buy_price pro home-01 (spot, buy_margin_percent=0)."""
|
||
dist = HOME01_DIST_VT_KWH if is_vt else HOME01_DIST_NT_KWH
|
||
if HOME01_BUY_MARGIN_PCT != 0.0:
|
||
if raw_ote >= 0:
|
||
energy = raw_ote * (1.0 + HOME01_BUY_MARGIN_PCT / 100.0)
|
||
else:
|
||
energy = raw_ote * (1.0 - HOME01_BUY_MARGIN_PCT / 100.0)
|
||
else:
|
||
energy = raw_ote
|
||
subtotal = (
|
||
energy
|
||
+ dist
|
||
+ HOME01_SYSTEM_SERVICES_KWH
|
||
+ HOME01_OTE_FEE_KWH
|
||
+ HOME01_BUY_MARGIN_FIXED_KWH
|
||
)
|
||
return subtotal * HOME01_VAT_MULTIPLIER
|
||
|
||
|
||
def build_buy_prices_96(raw_ote_96: Sequence[float], cfg: BuyPricingConfig) -> list[float]:
|
||
fixed_fees_kwh = cfg.distribution_kwh + cfg.other_fees_kwh
|
||
if cfg.mode == "spot_add_fixed":
|
||
if cfg.spot_add_fixed_kwh is None:
|
||
raise ValueError("Pro mode=spot_add_fixed chybí spot_add_fixed_kwh")
|
||
return [
|
||
(effective_buy_spot_add_fixed_kc_kwh(px, cfg.spot_add_fixed_kwh) + fixed_fees_kwh)
|
||
* cfg.vat_multiplier
|
||
for px in raw_ote_96
|
||
]
|
||
if cfg.mode == "spot_asym_pct":
|
||
if cfg.spot_asym_pct is None:
|
||
raise ValueError("Pro mode=spot_asym_pct chybí spot_asym_pct")
|
||
return [
|
||
(effective_buy_spot_asym_pct_kc_kwh(px, cfg.spot_asym_pct) + fixed_fees_kwh)
|
||
* cfg.vat_multiplier
|
||
for px in raw_ote_96
|
||
]
|
||
if cfg.mode == "nt_vt":
|
||
if cfg.nt_kwh is None or cfg.vt_kwh is None:
|
||
raise ValueError("Pro mode=nt_vt chybí NT/VT cena")
|
||
return [
|
||
(base + fixed_fees_kwh) * cfg.vat_multiplier
|
||
for base in buy_prices_96_nt_vt(
|
||
cfg.nt_kwh,
|
||
cfg.vt_kwh,
|
||
cfg.nt_from_hour,
|
||
cfg.nt_to_hour,
|
||
)
|
||
]
|
||
if cfg.mode == "home01":
|
||
return [
|
||
effective_buy_home01_kc_kwh(px, slot_is_vt_home01(t))
|
||
for t, px in enumerate(raw_ote_96)
|
||
]
|
||
if cfg.mode != "flat":
|
||
raise ValueError(f"Neznámý buy mode: {cfg.mode}")
|
||
return [(cfg.flat_kwh + fixed_fees_kwh) * cfg.vat_multiplier] * SLOTS_PER_DAY
|
||
|
||
|
||
def load_env_file(path: Path) -> None:
|
||
if not path.is_file():
|
||
return
|
||
for line in path.read_text(encoding="utf-8", errors="replace").splitlines():
|
||
line = line.strip()
|
||
if not line or line.startswith("#") or "=" not in line:
|
||
continue
|
||
k, _, v = line.partition("=")
|
||
k, v = k.strip(), v.strip().strip('"').strip("'")
|
||
if not k:
|
||
continue
|
||
if k not in os.environ or os.environ.get(k, "") == "":
|
||
os.environ[k] = v
|
||
|
||
|
||
def apply_auto_env_files() -> None:
|
||
"""Na produkci: /opt/ems-deploy/.env (když tam je docker-compose), pak app/.env nebo kořen repa."""
|
||
script = Path(__file__).resolve()
|
||
if len(script.parents) >= 3:
|
||
deploy_root = script.parents[2]
|
||
if (deploy_root / "docker-compose.yml").is_file():
|
||
load_env_file(deploy_root / ".env")
|
||
if len(script.parents) >= 2:
|
||
load_env_file(script.parents[1] / ".env")
|
||
|
||
|
||
def sync_pg_env_from_db_vars() -> None:
|
||
if not os.environ.get("PGUSER") and os.environ.get("DB_USER"):
|
||
os.environ["PGUSER"] = os.environ["DB_USER"]
|
||
if not os.environ.get("PGPASSWORD") and os.environ.get("DB_PASSWORD"):
|
||
os.environ["PGPASSWORD"] = os.environ["DB_PASSWORD"]
|
||
if not os.environ.get("PGDATABASE"):
|
||
os.environ["PGDATABASE"] = "ems"
|
||
|
||
|
||
def load_prices_csv(path: str) -> list[tuple[datetime, float]]:
|
||
out: list[tuple[datetime, float]] = []
|
||
with open(path, newline="", encoding="utf-8") as f:
|
||
first_line = f.readline()
|
||
f.seek(0)
|
||
if "interval_start" in first_line and "sell_raw_price_czk_kwh" in first_line:
|
||
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))
|
||
else:
|
||
from zoneinfo import ZoneInfo
|
||
|
||
prg = ZoneInfo("Europe/Prague")
|
||
r = csv.reader(f)
|
||
for row in r:
|
||
if len(row) < 3:
|
||
continue
|
||
date_s = row[0].strip()
|
||
time_s = row[1].strip()
|
||
price_s = row[2].strip()
|
||
if not date_s or not time_s or not price_s:
|
||
continue
|
||
ts = datetime.fromisoformat(f"{date_s}T{time_s}").replace(tzinfo=prg)
|
||
px = float(price_s)
|
||
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)
|
||
|
||
dsn = (
|
||
os.environ.get("DATABASE_URL")
|
||
or os.environ.get("EMS_DATABASE_URL")
|
||
or os.environ.get("POSTGRES_URL")
|
||
)
|
||
if dsn:
|
||
if dsn.startswith("postgres://"):
|
||
dsn = "postgresql://" + dsn[len("postgres://") :]
|
||
conn = psycopg2.connect(dsn)
|
||
else:
|
||
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: Sequence[float],
|
||
e_usable_wh: float,
|
||
p_batt_w: float,
|
||
site: SiteLimits,
|
||
soc_start_wh: float,
|
||
solver_time_limit_sec: 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)
|
||
batt_is_charging = pulp.LpVariable.dicts(
|
||
"batt_is_charging", range(SLOTS_PER_DAY), lowBound=0, upBound=1, cat="Binary"
|
||
)
|
||
grid_is_import = pulp.LpVariable.dicts(
|
||
"grid_is_import", range(SLOTS_PER_DAY), lowBound=0, upBound=1, cat="Binary"
|
||
)
|
||
|
||
prob += soc[0] == soc_start_wh
|
||
|
||
obj = []
|
||
for t in range(SLOTS_PER_DAY):
|
||
# Anti-loop constraints: a slot cannot both charge and discharge the battery,
|
||
# and it cannot both import from and export to the grid.
|
||
prob += ch[t] <= max_ch * batt_is_charging[t]
|
||
prob += dis[t] <= max_dis * (1 - batt_is_charging[t])
|
||
prob += gimp[t] <= max_imp * grid_is_import[t]
|
||
prob += gexp[t] <= max_exp * (1 - grid_is_import[t])
|
||
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[t] * gimp[t] / 1000.0)
|
||
|
||
prob += pulp.lpSum(obj)
|
||
|
||
solver_kwargs: dict[str, object] = {"msg": False}
|
||
if solver_time_limit_sec > 0:
|
||
solver_kwargs["timeLimit"] = solver_time_limit_sec
|
||
solver = pulp.PULP_CBC_CMD(**solver_kwargs)
|
||
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: Sequence[date],
|
||
px_day: dict[date, list[float]],
|
||
usable_kwh: float,
|
||
site: SiteLimits,
|
||
sell_margin_fixed: float,
|
||
sell_margin_pct: float,
|
||
buy_cfg: BuyPricingConfig,
|
||
summer_kwh: float,
|
||
winter_kwh: float,
|
||
load_kw: float,
|
||
shape: Sequence[float],
|
||
monthly_ed_kwh: Mapping[int, float] | None,
|
||
solver_time_limit_sec: float,
|
||
progress_every_days: int,
|
||
) -> 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
|
||
run_days = [d for d in days if d in px_day]
|
||
total_days = len(run_days)
|
||
started = perf_counter()
|
||
n_days = 0
|
||
if progress_every_days > 0:
|
||
limit_msg = (
|
||
f"{solver_time_limit_sec:g} s/den"
|
||
if solver_time_limit_sec > 0
|
||
else "bez limitu / den"
|
||
)
|
||
print(
|
||
f"[{usable_kwh:.1f} kWh] start: {total_days} dnů, CBC limit {limit_msg}",
|
||
flush=True,
|
||
)
|
||
for idx, d in enumerate(run_days, start=1):
|
||
if progress_every_days > 0 and (
|
||
idx == 1 or idx % progress_every_days == 0 or idx == total_days
|
||
):
|
||
elapsed_sec = perf_counter() - started
|
||
print(
|
||
f"[{usable_kwh:.1f} kWh] den {idx}/{total_days}: {d.isoformat()} "
|
||
f"(elapsed {elapsed_sec:.1f} s)",
|
||
flush=True,
|
||
)
|
||
raw = px_day[d]
|
||
p_sell = [effective_sell_kc_kwh(x, sell_margin_fixed, sell_margin_pct) for x in raw]
|
||
p_buy = build_buy_prices_96(raw, buy_cfg)
|
||
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,
|
||
p_buy,
|
||
e_wh,
|
||
p_batt,
|
||
site,
|
||
soc_state,
|
||
solver_time_limit_sec,
|
||
)
|
||
cash_total += cash
|
||
curt_total += curt
|
||
dis_total += dis
|
||
n_days += 1
|
||
if progress_every_days > 0:
|
||
print(
|
||
f"[{usable_kwh:.1f} kWh] hotovo za {perf_counter() - started:.1f} s",
|
||
flush=True,
|
||
)
|
||
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_* / DATABASE_URL)")
|
||
ap.add_argument(
|
||
"--no-auto-env",
|
||
action="store_true",
|
||
help="Nenačítej automaticky .env z ems-deploy/ ani app/",
|
||
)
|
||
ap.add_argument("--env-file", type=str, default="", help="Dodatečný soubor .env (po auto-env)")
|
||
ap.add_argument(
|
||
"--pg-host",
|
||
type=str,
|
||
default="",
|
||
help="Přepíše PGHOST (např. stejná IP jako EMS_DB_BIND ve compose)",
|
||
)
|
||
ap.add_argument(
|
||
"--price-csv",
|
||
type=str,
|
||
default="",
|
||
help="CSV buď s hlavičkou interval_start,sell_raw_price_czk_kwh, nebo legacy bez hlavičky: date,time,price",
|
||
)
|
||
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="Flat základní nákup Kč/kWh (když není spotový ani NT/VT režim)",
|
||
)
|
||
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(
|
||
"--buy-spot-add-fixed-kwh",
|
||
type=float,
|
||
default=None,
|
||
help="Základ nákupu = raw OTE + tento fixní adder Kč/kWh; pak se přičtou distribuce a ostatní poplatky",
|
||
)
|
||
ap.add_argument(
|
||
"--buy-spot-asym-pct",
|
||
type=float,
|
||
default=None,
|
||
help="Základ nákupu = raw OTE × (1 + p/100) pro raw >= 0, raw OTE × (1 - p/100) pro raw < 0",
|
||
)
|
||
ap.add_argument(
|
||
"--buy-home-01",
|
||
action="store_true",
|
||
help="Nákup jako home-01 z živé DB (spot ×1.09/×0.91, dist NT/VT HDO, SS+OTE, DPH 21 %)",
|
||
)
|
||
ap.add_argument(
|
||
"--buy-distribution-kwh",
|
||
type=float,
|
||
default=0.0,
|
||
help="Fixně přičtená distribuční složka Kč/kWh ke každému nákupnímu slotu",
|
||
)
|
||
ap.add_argument(
|
||
"--buy-other-fees-kwh",
|
||
type=float,
|
||
default=0.0,
|
||
help="Fixně přičtené ostatní poplatky Kč/kWh (OTE, systémové služby apod.) ke každému nákupnímu slotu",
|
||
)
|
||
ap.add_argument(
|
||
"--buy-vat-multiplier",
|
||
type=float,
|
||
default=1.0,
|
||
help="Násobitel DPH aplikovaný na finální nákupní cenu po přičtení distribuce a ostatních poplatků (např. 1.21)",
|
||
)
|
||
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)
|
||
ap.add_argument("--c-rate", type=float, default=0.5)
|
||
ap.add_argument(
|
||
"--solver-time-limit-sec",
|
||
type=float,
|
||
default=60.0,
|
||
help="CBC time limit na jeden den; 0 = bez limitu",
|
||
)
|
||
ap.add_argument(
|
||
"--progress-every-days",
|
||
type=int,
|
||
default=1,
|
||
help="Po kolika dnech vytisknout průběh; 0 = tichý režim",
|
||
)
|
||
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)
|
||
base_buy_modes = [
|
||
args.buy_nt_kwh is not None,
|
||
args.buy_spot_add_fixed_kwh is not None,
|
||
args.buy_spot_asym_pct is not None,
|
||
args.buy_home_01,
|
||
]
|
||
if sum(base_buy_modes) > 1:
|
||
ap.error(
|
||
"Zvol jen jeden režim nákupu: flat, NT/VT, --buy-spot-add-fixed-kwh, "
|
||
"--buy-spot-asym-pct nebo --buy-home-01"
|
||
)
|
||
if args.buy_vat_multiplier <= 0:
|
||
ap.error("--buy-vat-multiplier musí být > 0")
|
||
if args.solver_time_limit_sec < 0:
|
||
ap.error("--solver-time-limit-sec musí být >= 0")
|
||
if args.progress_every_days < 0:
|
||
ap.error("--progress-every-days musí být >= 0")
|
||
for hour_arg, hour_value in (("nt-from-hour", args.nt_from_hour), ("nt-to-hour", args.nt_to_hour)):
|
||
if not (0 <= hour_value <= 23):
|
||
ap.error(f"--{hour_arg} musí být v rozsahu 0..23")
|
||
if args.db:
|
||
if not args.no_auto_env:
|
||
apply_auto_env_files()
|
||
if args.env_file:
|
||
load_env_file(Path(args.env_file))
|
||
sync_pg_env_from_db_vars()
|
||
if args.pg_host:
|
||
os.environ["PGHOST"] = args.pg_host
|
||
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,
|
||
)
|
||
|
||
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_home_01:
|
||
buy_cfg = BuyPricingConfig(mode="home01")
|
||
elif args.buy_spot_add_fixed_kwh is not None:
|
||
buy_cfg = BuyPricingConfig(
|
||
mode="spot_add_fixed",
|
||
spot_add_fixed_kwh=args.buy_spot_add_fixed_kwh,
|
||
distribution_kwh=args.buy_distribution_kwh,
|
||
other_fees_kwh=args.buy_other_fees_kwh,
|
||
vat_multiplier=args.buy_vat_multiplier,
|
||
)
|
||
elif args.buy_spot_asym_pct is not None:
|
||
buy_cfg = BuyPricingConfig(
|
||
mode="spot_asym_pct",
|
||
spot_asym_pct=args.buy_spot_asym_pct,
|
||
distribution_kwh=args.buy_distribution_kwh,
|
||
other_fees_kwh=args.buy_other_fees_kwh,
|
||
vat_multiplier=args.buy_vat_multiplier,
|
||
)
|
||
elif args.buy_nt_kwh is not None:
|
||
vt = args.buy_nt_kwh + args.buy_vt_surcharge_kwh
|
||
buy_cfg = BuyPricingConfig(
|
||
mode="nt_vt",
|
||
nt_kwh=args.buy_nt_kwh,
|
||
vt_kwh=vt,
|
||
nt_from_hour=args.nt_from_hour,
|
||
nt_to_hour=args.nt_to_hour,
|
||
distribution_kwh=args.buy_distribution_kwh,
|
||
other_fees_kwh=args.buy_other_fees_kwh,
|
||
vat_multiplier=args.buy_vat_multiplier,
|
||
)
|
||
else:
|
||
buy_cfg = BuyPricingConfig(
|
||
mode="flat",
|
||
flat_kwh=args.buy_vat_kwh,
|
||
distribution_kwh=args.buy_distribution_kwh,
|
||
other_fees_kwh=args.buy_other_fees_kwh,
|
||
vat_multiplier=args.buy_vat_multiplier,
|
||
)
|
||
|
||
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,
|
||
buy_cfg,
|
||
args.pv_daily_kwh_summer,
|
||
args.pv_daily_kwh_winter,
|
||
args.load_kw,
|
||
shape,
|
||
monthly_ed,
|
||
args.solver_time_limit_sec,
|
||
args.progress_every_days,
|
||
)
|
||
results.append((kwh, r))
|
||
|
||
baseline_kwh = min(args.battery_kwh)
|
||
base = dict(results)[baseline_kwh]
|
||
|
||
print("Parametry: prodej = OTE + sell_margin_fixed (+ %)")
|
||
if buy_cfg.mode == "nt_vt":
|
||
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í)"
|
||
)
|
||
elif buy_cfg.mode == "spot_add_fixed":
|
||
print(f" Nákup = raw OTE + {args.buy_spot_add_fixed_kwh} Kč/kWh")
|
||
elif buy_cfg.mode == "spot_asym_pct":
|
||
print(
|
||
f" Nákup = raw OTE × (1 + {args.buy_spot_asym_pct}/100) pro raw >= 0, "
|
||
f"raw OTE × (1 - {args.buy_spot_asym_pct}/100) pro raw < 0"
|
||
)
|
||
elif buy_cfg.mode == "home01":
|
||
print(
|
||
" Nákup = home-01 (MCP): raw OTE ×(1+9%) / ×(1-9%) + dist NT/VT dle HDO + "
|
||
f"{HOME01_SYSTEM_SERVICES_KWH} SS + {HOME01_OTE_FEE_KWH} OTE, ×{HOME01_VAT_MULTIPLIER} DPH"
|
||
)
|
||
print(
|
||
f" distribuce NT {HOME01_DIST_NT_KWH} / VT {HOME01_DIST_VT_KWH} Kč/kWh; "
|
||
"VT 09–10, 12–13, 16–17, 20–21; prodej typicky sell_margin_fixed -0.30"
|
||
)
|
||
else:
|
||
print(f" Nákup = flat {args.buy_vat_kwh} Kč/kWh")
|
||
if args.buy_distribution_kwh or args.buy_other_fees_kwh:
|
||
print(
|
||
f" Fixní add-on k nákupu: distribuce {args.buy_distribution_kwh} Kč/kWh, "
|
||
f"ostatní poplatky {args.buy_other_fees_kwh} Kč/kWh"
|
||
)
|
||
if args.buy_vat_multiplier != 1.0:
|
||
print(f" DPH násobitel na finální nákupní cenu: {args.buy_vat_multiplier}")
|
||
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)")
|
||
limit_msg = (
|
||
f"{args.solver_time_limit_sec:g} s/den"
|
||
if args.solver_time_limit_sec > 0
|
||
else "bez limitu / den"
|
||
)
|
||
print(f" Solver: CBC, limit {limit_msg}, progress every {args.progress_every_days} dnů")
|
||
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()
|