vpp-ai-platform/skills-py/vpp_skills/forecast.py

"""Forecast skills v1: load / PV / price with quantile intervals (docs/05 §2.1).

Baseline method, deliberately simple and fully deterministic:

* point forecast (P50): exponentially-weighted mean of recent history days of
  the same day-type (weekday vs weekend), per 15-min interval;
* interval (P10/P90): empirical quantiles of *backtest residuals* — the same
  method is re-run on each history day using only the days before it, and the
  pooled (actual − predicted) residuals give the band. This is split-conformal
  calibration in its simplest form, so coverage on held-out days tracks the
  nominal 80% instead of being an optimistic guess (docs/12 L2: 区间校准单独考核).

Better models (statsforecast, gradient boosting) plug in behind the same
contract; the eval harness (packages/evals) is what decides whether they win.
Exogenous inputs are accepted by the contract but unused by v1 — recorded here
so nobody mistakes v1 for weather-aware.
"""

from __future__ import annotations

from datetime import date, datetime, timezone
from typing import Callable

import numpy as np

from vpp_contracts.forecast_bundle import ForecastBundle
from vpp_contracts.forecast_request import ForecastRequest

from . import SKILL_VERSIONS
from .numeric import INTERVALS, SCALE_MW, SCALE_PRICE, curve, to_array

MODEL_NAME = {"LOAD": "load-forecast", "PV": "pv-forecast", "PRICE": "price-forecast"}
UNIT_FOR_KIND = {"LOAD": "mw", "PV": "mw", "PRICE": "yuan_per_mwh"}

LOOKBACK_DAYS = 14
HALF_LIFE_DAYS = 3.0
MIN_BACKTEST_HISTORY = 3
NOMINAL_LOWER, NOMINAL_UPPER = 0.10, 0.90


def _is_weekend(d: str) -> bool:
    return date.fromisoformat(d).weekday() >= 5


def _point_forecast(days: list[tuple[str, np.ndarray]], target: str) -> np.ndarray:
    """EWM over the most recent LOOKBACK_DAYS same-day-type days (fallback: all days)."""
    same = [(d, v) for d, v in days if _is_weekend(d) == _is_weekend(target)]
    pool = (same if len(same) >= 2 else days)[-LOOKBACK_DAYS:]
    ages = np.array([(date.fromisoformat(target) - date.fromisoformat(d)).days for d, _ in pool], float)
    w = np.power(0.5, ages / HALF_LIFE_DAYS)
    stack = np.stack([v for _, v in pool])
    return (w[:, None] * stack).sum(axis=0) / w.sum()


def _residual_quantiles(days: list[tuple[str, np.ndarray]], relative: bool) -> tuple[float, float]:
    """Pooled backtest residuals. Relative (actual/pred − 1) for quantities whose
    spread scales with level (load, PV — a night-time zero must not shrink the
    noon band); additive for prices, which can sit near or below zero."""
    pairs: list[tuple[np.ndarray, np.ndarray]] = []
    for i in range(MIN_BACKTEST_HISTORY, len(days)):
        d, actual = days[i]
        pairs.append((actual, _point_forecast(days[:i], d)))
    if not pairs:
        # Too little history to backtest: fall back to day-to-day differences.
        pairs = [(days[i][1], days[i - 1][1]) for i in range(1, len(days))]
    actual = np.concatenate([a for a, _ in pairs])
    pred = np.concatenate([p for _, p in pairs])
    if relative:
        floor = 0.05 * max(float(pred.max()), 1e-9)
        mask = pred > floor
        pooled = actual[mask] / pred[mask] - 1.0 if mask.any() else np.zeros(1)
    else:
        pooled = actual - pred
    lo, hi = np.quantile(pooled, [NOMINAL_LOWER, NOMINAL_UPPER])
    return float(lo), float(hi)


def forecast(
    req: ForecastRequest,
    clock: Callable[[], datetime] = lambda: datetime.now(timezone.utc),
) -> ForecastBundle:
    kind = str(req.kind.value)
    if req.unit.value != UNIT_FOR_KIND[kind]:
        raise ValueError(f"{kind} forecast unit must be {UNIT_FOR_KIND[kind]}, got {req.unit.value}")

    days = [(c.date, to_array(v.root for v in c.values)) for c in req.history]
    dates = [d for d, _ in days]
    if dates != sorted(dates) or len(set(dates)) != len(dates):
        raise ValueError("history must be ascending by date with no duplicates")
    if dates[-1] >= req.market_date:
        raise ValueError("history must be strictly before market_date")

    p50 = _point_forecast(days, req.market_date)
    relative = kind in ("LOAD", "PV")
    lo, hi = _residual_quantiles(days, relative)
    if relative:
        p10, p90 = p50 * (1.0 + lo), p50 * (1.0 + hi)
    else:
        p10, p90 = p50 + lo, p50 + hi

    if kind in ("LOAD", "PV"):
        p10, p50, p90 = (np.clip(x, 0.0, None) for x in (p10, p50, p90))
    if kind == "PV":
        # Intervals that were zero on every history day (night) stay exactly zero.
        night = np.all(np.stack([v for _, v in days]) == 0.0, axis=0)
        for x in (p10, p50, p90):
            x[night] = 0.0

    scale = SCALE_PRICE if kind == "PRICE" else SCALE_MW
    assert p50.shape == (INTERVALS,)
    bundle = {
        "id": f"fc-{kind.lower()}-{req.market_date}",
        "kind": kind,
        "market_date": req.market_date,
        "unit": req.unit.value,
        "quantiles": {
            "p10": curve(p10, req.market_date, scale),
            "p50": curve(p50, req.market_date, scale),
            "p90": curve(p90, req.market_date, scale),
        },
        "model": {"name": MODEL_NAME[kind], "version": SKILL_VERSIONS[MODEL_NAME[kind]]},
        "features_snapshot_ref": req.features_snapshot_ref,
        "generated_at": clock().isoformat().replace("+00:00", "Z"),
    }
    return ForecastBundle.model_validate(bundle)