ems/scripts/harness/hu1_bess_study.py

#!/usr/bin/env python3
"""
Studie HU1 (Hulín BESS): 128 kWh / 36 kW baterie, 2×20 kW Deye (AC 40 kW).

Ekonomika čistého BESS (bez PV a spotřeby) nad reálným obdobím BA81 auditu:
  - nákup FIXNÍ jako BA81 (efektivní buy site 3 — dle záměru smlouvy),
  - prodej SPOT (efektivní sell site 3 jako proxy; HU1 marže dle budoucí smlouvy),
  - limity: baterie 36 kW (BMS/střídač), AC 40 kW, import 43 kW, export 42 kW,
    block_export_on_negative_sell = true (konfigurace site 5 v DB),
  - volitelný EDC sdílecí kanál z BA81: v slotech, kdy má BA81 přebytek při
    sell<0, lze nabíjet za pouhou distribuci (parametr; citlivost 1.0/1.5/2.0
    Kč/kWh). Množství = skutečný přebytek BA81 (pv − load) z auditu.

Výstup: Kč/den bez sdílení a se sdílením — horní mez (perfect hindsight).
POZOR: EDC sdílení zatím v EMS NENÍ implementováno — viz závěr studie.
"""

from __future__ import annotations

import asyncio
import os
import sys
from dataclasses import dataclass

import asyncpg
import pulp

INTERVAL_H = 0.25
WINDOW_DAYS = 7

BAT_USABLE_WH = 128_000.0
BAT_POWER_W = 36_000.0
AC_W = 40_000.0
IMP_W = 43_000.0
EXP_W = 42_000.0
MIN_PCT, MAX_PCT = 10.0, 100.0
EFF = 0.95
DEG_CZK_KWH = 0.30
BLOCK_NEG = True
DIST_SENSITIVITY = [None, 2.0, 1.5, 1.0]  # None = bez sdílení


@dataclass
class Slot:
    buy: float
    sell: float
    share_wh: float  # sdílitelný přebytek BA81 (jen sloty sell<0, jinak 0)


def _dsn() -> str:
    return os.environ.get(
        "EMS_DB_DSN", "postgresql://ems_user:dev_password@10.200.200.1:5432/ems"
    )


async def _load(conn: asyncpg.Connection, price_site: int = 3) -> list[Slot]:
    rows = await conn.fetch(
        """
        select case when $1 = 5 then p2.effective_buy_price_czk_kwh
                    else p.effective_buy_price_czk_kwh end as buy,
               case when $1 = 5 then p2.effective_sell_price_czk_kwh
                    else p.effective_sell_price_czk_kwh end as sell,
               -- POTENCIÁL: BA81 při sell<0 škrtí (81 % výroby v datech chybí)
               case when p.effective_sell_price_czk_kwh < 0
                 then greatest(0,
                        greatest(coalesce(a.actual_pv_production_wh,0),
                                 coalesce(fc.fc_a_wh,0) + coalesce(fc.fc_b_wh,0))
                        - coalesce(a.actual_load_consumption_wh,0))
                 else 0 end as share_wh
        from ems.audit_interval a
        join ems.vw_site_effective_price p
          on p.site_id = a.site_id and p.interval_start = a.interval_start
        left join ems.vw_site_effective_price p2
          on p2.site_id = 5 and p2.interval_start = a.interval_start
        left join lateral (
          select
            sum(power_w) filter (where pa.controllable) * 0.25 as fc_a_wh,
            sum(power_w) filter (where not pa.controllable) * 0.25 as fc_b_wh
          from (
            select distinct on (fpr.pv_array_id) fpi2.power_w, fpr.pv_array_id
            from ems.forecast_pv_interval fpi2
            join ems.forecast_pv_run fpr on fpr.id = fpi2.run_id
            where fpi2.interval_start = a.interval_start
            order by fpr.pv_array_id, fpr.created_at desc
          ) x
          join ems.asset_pv_array pa on pa.id = x.pv_array_id and pa.site_id = a.site_id
        ) fc on true
        where a.site_id = 3 and a.actual_load_consumption_wh is not null
          and p2.effective_buy_price_czk_kwh is not null
        order by a.interval_start
        """,
        price_site,
    )
    return [Slot(float(r["buy"]), float(r["sell"]), float(r["share_wh"])) for r in rows]


def solve_window(slots: list[Slot], soc0: float, dist: float | None) -> tuple[float, float]:
    n = len(slots)
    prob = pulp.LpProblem("hu1", pulp.LpMinimize)
    mc = min(BAT_POWER_W, AC_W) * INTERVAL_H
    mi = IMP_W * INTERVAL_H
    me = min(EXP_W, AC_W) * INTERVAL_H
    smin = MIN_PCT / 100 * BAT_USABLE_WH
    smax = MAX_PCT / 100 * BAT_USABLE_WH

    gi = [pulp.LpVariable(f"gi{t}", 0, mi) for t in range(n)]
    ge = [pulp.LpVariable(f"ge{t}", 0, me) for t in range(n)]
    bc = [pulp.LpVariable(f"bc{t}", 0, mc) for t in range(n)]
    bd = [pulp.LpVariable(f"bd{t}", 0, mc) for t in range(n)]
    sh = [
        pulp.LpVariable(f"sh{t}", 0, slots[t].share_wh if dist is not None else 0)
        for t in range(n)
    ]
    soc = [pulp.LpVariable(f"s{t}", smin, smax) for t in range(n)]
    y = [pulp.LpVariable(f"y{t}", cat="Binary") for t in range(n)]

    for t in range(n):
        s = slots[t]
        # BESS bez lokální spotřeby/PV: nabíjení ze sítě/sdílení, vybíjení do exportu
        prob += gi[t] + sh[t] + bd[t] == bc[t] + ge[t]
        prev = soc0 if t == 0 else soc[t - 1]
        prob += soc[t] == prev + bc[t] * EFF - bd[t] / EFF
        prob += gi[t] + sh[t] <= mi * y[t]
        prob += ge[t] <= me * (1 - y[t])
        if s.sell < 0 and BLOCK_NEG:
            prob += ge[t] == 0
        if s.buy < 0:
            prob += ge[t] == 0

    avg_buy = sum(s.buy for s in slots) / n
    cash = pulp.lpSum(
        gi[t] / 1000 * slots[t].buy
        + (sh[t] / 1000 * dist if dist is not None else 0)
        - ge[t] / 1000 * slots[t].sell
        for t in range(n)
    )
    deg = pulp.lpSum(0.5 * (bc[t] + bd[t]) / 1000 * DEG_CZK_KWH for t in range(n))
    prob += cash + deg - soc[n - 1] / 1000 * max(0.0, avg_buy)

    solver = pulp.HiGHS_CMD(msg=False, timeLimit=30) if pulp.HiGHS_CMD().available() else pulp.PULP_CBC_CMD(msg=False)
    prob.solve(solver)
    if pulp.LpStatus[prob.status] != "Optimal":
        raise RuntimeError(pulp.LpStatus[prob.status])
    return float(pulp.value(cash)) + float(pulp.value(deg)), float(soc[n - 1].value())


async def main() -> None:
    conn = await asyncpg.connect(_dsn())
    try:
        import os as _os
        price_site = int(_os.environ.get('HU1_PRICE_SITE', '3'))
        slots = await _load(conn, price_site)
    finally:
        await conn.close()
    days_total = len(slots) // 96
    share_total = sum(s.share_wh for s in slots) / 1000
    print(f"# HU1 BESS studie — 128 kWh / 36 kW / AC 40 kW; ceny site {price_site} ({'fixní nákup BA81' if price_site==3 else 'SPOT nákup i prodej (site 5)'})")
    print(f"# Období: {days_total} dní (BA81 audit); sdílitelný přebytek BA81 při sell<0: {share_total:.0f} kWh\n")
    for dist in DIST_SENSITIVITY:
        total = 0.0
        soc = 0.3 * BAT_USABLE_WH
        days = 0
        i = 0
        while i + 96 <= len(slots):
            win = slots[i : i + WINDOW_DAYS * 96]
            if len(win) < 96:
                break
            cost, soc = solve_window(win, soc, dist)
            total += cost
            days += len(win) // 96
            i += len(win)
        per_day = -total / max(1, days)
        label = "bez EDC sdílení" if dist is None else f"EDC sdílení, distribuce {dist:.1f} Kč/kWh"
        print(f"  {label:<42} výnos {-total:>8.0f} Kč  = {per_day:>7.2f} Kč/den  (rok ~{per_day*365*0.6:.0f}–{per_day*365:.0f} Kč)")
    print("\nPozn.: horní mez (perfect hindsight); jaro = nejsilnější sezóna pro spot spready.")


if __name__ == "__main__":
    sys.exit(asyncio.run(main()))