Day-ahead PV forecast for Germany

import os
import sys

# notebooks/ is the documented home; run from the project root either way
if os.path.basename(os.getcwd()) == "notebooks":
    os.chdir("..")
sys.path.insert(0, os.getcwd())

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.inspection import permutation_importance
from sklearn.metrics import mean_absolute_error, mean_pinball_loss

from src.build_baselines import all_baselines
from src.build_features import HORIZON_H, feature_columns
from src.evaluate import CAPACITY_MW, score_table
from src.fetch import SITES

SEED = 0            # fixed everywhere so the README numbers reproduce
TRAIN, TEST = "2023", "2024"

# rows the lag features cannot supply; never dropna on the feature columns
# (k_cs_lag24 is NaN at night by design and HistGradientBoosting handles NaN)
NEED = ["pv_actual", "pv_lag24", "pv_lag168", "pv_roll7d"]

pd.set_option("display.width", 120, "display.max_columns", 30)
plt.rcParams.update({"figure.dpi": 110, "axes.grid": True, "grid.alpha": 0.3})

feat = pd.read_parquet("data/cached/features.parquet")

print(f"{feat.shape[0]:,} hours x {feat.shape[1]} columns")
print(f"{feat.index.min()}  ->  {feat.index.max()}  ({feat.index.tz})")
print(f"horizon {HORIZON_H} h, nMAE normalised by {CAPACITY_MW:,} MW")

cols = feature_columns(feat)
FORBIDDEN = {"pv_actual", "pv_forecast", "k_cs", "is_day", "load"}
assert not (FORBIDDEN & set(cols)), FORBIDDEN & set(cols)

per_city = [c for c in cols if any(c.endswith(f"_{s}") for s in SITES)]
national = [c for c in cols if c not in per_city]

print(f"{len(cols)} features = {len(national)} national/calendar "
      f"+ {len(per_city)} per-city ({len(per_city) // len(SITES)} x {len(SITES)} cities)")
print(f"\nnational and calendar:\n  {national}")
print(f"\nper city ({', '.join(SITES)}):")
print("  " + str(sorted({c.rsplit('_', 1)[0] for c in per_city})))
print(f"\nheld out of the model: {sorted(FORBIDDEN)}")

train = feat.loc[TRAIN].dropna(subset=NEED)
test = feat.loc[TEST].dropna(subset=NEED)

for name, part in [("train", train), ("test", test)]:
    print(f"{name:5s} {len(part):5,} h   {part.index.min().date()} -> "
          f"{part.index.max().date()}   {part['is_day'].sum():,} daylight h")
assert train.index.max() < test.index.min(), "train must end before test starts"

preds = all_baselines(feat).loc[test.index]
print(score_table(test, preds))

gbr = HistGradientBoostingRegressor(random_state=SEED)
gbr.fit(train[cols], train["pv_actual"])

preds["model"] = gbr.predict(test[cols]).clip(0)
print(gbr)

# overfit check: train and test error side by side
for name, part in [("train", train), ("test", test)]:
    p = pd.Series(gbr.predict(part[cols]), index=part.index).clip(0)
    day = part["is_day"]
    print(f"{name:5s} MAE {mean_absolute_error(part['pv_actual'], p):6.0f} MW   "
          f"MAE_day {mean_absolute_error(part['pv_actual'][day], p[day]):6.0f} MW")

tab = score_table(test, preds)

cs = tab.loc["clearsky_persistence", "MAE_day"]
tso = tab.loc["tso_forecast", "MAE_day"]
tab["skill_vs_cs"] = (1 - tab["MAE_day"] / cs).round(3)
tab["gap_to_tso_%"] = (100 * (tab["MAE_day"] / tso - 1)).round(1)

print(f"{TEST} test set, {tab.attrs['n_hours']:,} hours scored "
      f"(identical rows for every method)\n")
print(tab)

err_day = {}
for name in ["clearsky_persistence", "tso_forecast", "model"]:
    e = (preds[name] - test["pv_actual"])[test["is_day"]].dropna()
    err_day[name] = {"bias": e.mean(), "MAE": e.abs().mean(),
                     "dispersion (MAE of e - bias)": (e - e.mean()).abs().mean()}

print(f"daylight hours only, mean actual "
      f"{test['pv_actual'][test['is_day']].mean():.0f} MW\n")
print(pd.DataFrame(err_day).T.round(0))

peak = feat["pv_actual"].resample("MS").max()
roll = feat["pv_actual"].rolling(24 * 30, min_periods=24 * 7).max()

fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(roll.index, roll, lw=1, label="30-day rolling max of pv_actual")
ax.plot(peak.index, peak, "o-", ms=4, lw=1, alpha=0.7, label="monthly peak")
ax.axvline(pd.Timestamp(f"{TEST}-01-01", tz="UTC"), color="k", ls="--", lw=1)
ax.text(pd.Timestamp(f"{TEST}-01-05", tz="UTC"), 5000, " test period", fontsize=9)
ax.set(title="Installed capacity is a moving target",
       ylabel="MW", xlabel="")
ax.legend(loc="upper left")
plt.tight_layout()
plt.show()

s23 = feat.loc[TRAIN, "pv_actual"].max()
s24 = feat.loc[TEST, "pv_actual"].max()
print(f"peak generation {TRAIN}: {s23:,.0f} MW   {TEST}: {s24:,.0f} MW   "
      f"(+{100 * (s24 / s23 - 1):.0f}%)")

cap = feat["pv_actual"].shift(HORIZON_H).rolling(24 * 30, min_periods=24 * 14).max()
assert "cap" not in cols, "capacity proxy must not become a feature"

g = feat.assign(cap=cap).dropna(subset=NEED + ["cap"])
tr_r, te_r = g.loc[TRAIN], g.loc[TEST]

ratio = (tr_r["pv_actual"] / tr_r["cap"]).clip(0, 1.5)
gbr_r = HistGradientBoostingRegressor(random_state=SEED).fit(tr_r[cols], ratio)
model_scaled = (pd.Series(gbr_r.predict(te_r[cols]), index=te_r.index)
                * te_r["cap"]).clip(0)

cmp = preds.copy()
cmp["model_capacity_scaled"] = model_scaled
tab2 = score_table(test, cmp)
tab2["bias_day"] = [
    (cmp[n] - test["pv_actual"])[test["is_day"]].mean().round(0) for n in cmp.columns
]
print(tab2.loc[["clearsky_persistence", "model", "model_capacity_scaled",
                "tso_forecast"], ["MAE_day", "RMSE_day", "nMAE_day_%", "bias_day"]])

METHODS = ["clearsky_persistence", "model", "tso_forecast"]
COLOURS = {"clearsky_persistence": "tab:grey", "model": "tab:blue",
           "tso_forecast": "tab:orange"}

d = (pd.concat([test[["pv_actual", "is_day", "cloud_cover_de"]], preds], axis=1)
       .dropna(subset=list(preds.columns) + ["pv_actual"]))
for name in METHODS:
    d[f"e_{name}"] = (d[name] - d["pv_actual"]).abs()
err_cols = [f"e_{name}" for name in METHODS]

print(f"{len(d):,} rows, matching the {tab.attrs['n_hours']:,} scored above")


def style(ax, title, ylabel, xlabel="", legend=METHODS):
    """legend=None keeps whatever labels the artists already carry."""
    ax.set(title=title, ylabel=ylabel, xlabel=xlabel)
    ax.legend(legend) if legend else ax.legend()

day_rows = d[d["is_day"]]
summer = day_rows[(day_rows.index.month >= 4) & (day_rows.index.month <= 9)]
weekly = summer.groupby(pd.Grouper(freq="W-MON"))["cloud_cover_de"].agg(["mean", "count"])
weekly = weekly[weekly["count"] >= 90]

picks = [("clearest", weekly["mean"].idxmin(), weekly["mean"].min()),
         ("cloudiest", weekly["mean"].idxmax(), weekly["mean"].max())]

fig, axes = plt.subplots(2, 1, figsize=(11, 7), sharey=True)
for ax, (label, end, cc) in zip(axes, picks):
    w = d.loc[end - pd.Timedelta(days=7):end]
    ax.plot(w.index, w["pv_actual"], color="k", lw=1.8, label="actual")
    for name in METHODS:
        ax.plot(w.index, w[name], lw=1.2, alpha=0.85, color=COLOURS[name], label=name)
    ax.set(title=f"{label} week of Apr-Sep {TEST} "
                 f"({(end - pd.Timedelta(days=7)).date()} to {end.date()}, "
                 f"mean daylight cloud cover {cc:.0f}%)",
           ylabel="MW")
    ax.legend(ncol=4, fontsize=8)
plt.tight_layout()
plt.show()

by_hour = d.groupby(d.index.hour)[err_cols].mean()

# relative: daylight only, since dividing by a near-zero night mean explodes
mean_by_hour = day_rows.groupby(day_rows.index.hour)["pv_actual"].mean()
rel_hour = 100 * day_rows.groupby(day_rows.index.hour)[err_cols].mean().div(
    mean_by_hour, axis=0)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
by_hour.plot(ax=axes[0], color=[COLOURS[m] for m in METHODS])
style(axes[0], f"Absolute error by hour, {TEST}", "MAE [MW]", "hour (UTC)")
rel_hour.plot(ax=axes[1], color=[COLOURS[m] for m in METHODS])
style(axes[1], f"Relative error by hour (daylight), {TEST}",
      "MAE / mean output [%]", "hour (UTC)")
plt.tight_layout()
plt.show()

daily = d[err_cols].resample("D").mean()
daily_pv = d["pv_actual"].resample("D").mean()
smooth = daily.rolling(30, center=True, min_periods=10).mean()
rel_smooth = (100 * daily.div(daily_pv, axis=0)).rolling(
    30, center=True, min_periods=10).mean()

fig, axes = plt.subplots(2, 1, figsize=(11, 7))
smooth.plot(ax=axes[0], color=[COLOURS[m] for m in METHODS])
style(axes[0], f"Absolute error through {TEST} (30-day running mean)", "MAE [MW]")
rel_smooth.plot(ax=axes[1], color=[COLOURS[m] for m in METHODS])
style(axes[1], f"Relative error through {TEST} (30-day running mean)",
      "MAE / mean output [%]")
plt.tight_layout()
plt.show()

monthly = day_rows.groupby(day_rows.index.month)[err_cols].mean()
monthly.columns = METHODS
monthly.index.name = "month"
monthly["mean_output"] = day_rows.groupby(day_rows.index.month)["pv_actual"].mean()
monthly["model_rel_%"] = 100 * monthly["model"] / monthly["mean_output"]
print(f"daylight MAE by month, {TEST} [MW]\n")
print(monthly.round(0))

by_day = day_rows.resample("D")[err_cols + ["pv_actual", "cloud_cover_de"]].mean()
worst = by_day["e_model"].dropna().nlargest(8)
detail = by_day.loc[worst.index, ["e_model", "e_tso_forecast", "pv_actual",
                                  "cloud_cover_de"]]
detail.columns = ["model_MAE", "tso_MAE", "mean_actual", "cloud_cover_%"]
print(f"worst 8 days for the model, daylight MAE [MW]\n")
print(detail.round(0))

bad = worst.index[0]
w = d.loc[bad.strftime("%Y-%m-%d")]

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(w.index, w["pv_actual"], color="k", lw=2, label="actual")
for name in METHODS:
    ax.plot(w.index, w[name], lw=1.4, color=COLOURS[name], label=name)
ax.fill_between(w.index, w["model"], w["pv_actual"], alpha=0.15, color="tab:blue")
style(ax, f"Worst day for the model: {bad.date()} "
          f"(daylight MAE {worst.iloc[0]:.0f} MW)", "MW", "hour (UTC)",
      legend=["actual"] + METHODS)
plt.tight_layout()
plt.show()

def fit_quantiles(X_train, y_train, X_test, scale=None, taus=(0.1, 0.5, 0.9)):
    """One fit per quantile; `scale` multiplies predictions back to MW."""
    out = {}
    for tau in taus:
        m = HistGradientBoostingRegressor(loss="quantile", quantile=tau,
                                          random_state=SEED)
        m.fit(X_train, y_train)
        p = pd.Series(m.predict(X_test), index=X_test.index)
        out[f"p{int(tau * 100)}"] = (p if scale is None else p * scale).clip(0)
    q = pd.DataFrame(out)
    # separate fits can cross; enforce p10 <= p50 <= p90
    q["p50"] = q[["p10", "p50"]].max(axis=1)
    q["p90"] = q[["p50", "p90"]].max(axis=1)
    return q


def report(q):
    y = test["pv_actual"].reindex(q.index)
    day = test["is_day"].reindex(q.index)
    inside = (y >= q["p10"]) & (y <= q["p90"])
    return {
        "pinball_p10": mean_pinball_loss(y[day], q["p10"][day], alpha=0.1),
        "pinball_p50": mean_pinball_loss(y[day], q["p50"][day], alpha=0.5),
        "pinball_p90": mean_pinball_loss(y[day], q["p90"][day], alpha=0.9),
        "coverage_%": 100 * inside[day].mean(),
        "below_p10_%": 100 * (y < q["p10"])[day].mean(),
        "above_p90_%": 100 * (y > q["p90"])[day].mean(),
        "mean_width": (q["p90"] - q["p10"])[day].mean(),
    }


q_raw = fit_quantiles(train[cols], train["pv_actual"], test[cols])
q_scaled = fit_quantiles(tr_r[cols], ratio, te_r[cols], scale=te_r["cap"])

print("daylight hours; nominal coverage 80%, 10% in each tail\n")
print(pd.DataFrame({"raw target": report(q_raw),
                    "capacity-scaled": report(q_scaled)}).T.round(1))

end = picks[1][1]                       # the cloudy week, where the band matters
w = q_scaled.loc[end - pd.Timedelta(days=5):end]
y = test["pv_actual"].reindex(w.index)

fig, ax = plt.subplots(figsize=(11, 4))
ax.fill_between(w.index, w["p10"], w["p90"], alpha=0.3, color="tab:blue",
                label="P10-P90")
ax.plot(w.index, w["p50"], color="tab:blue", lw=1.4, label="P50")
ax.plot(y.index, y, color="k", lw=1.8, label="actual")
ax.plot(w.index, preds["tso_forecast"].reindex(w.index), color="tab:orange",
        lw=1.2, alpha=0.8, label="tso_forecast")
inside = ((y >= w["p10"]) & (y <= w["p90"]))
ax.set(title=f"Capacity-scaled quantiles, cloudy week of {TEST} "
             f"({100 * inside.mean():.0f}% of hours inside the band)",
       ylabel="MW", xlabel="")
ax.legend(ncol=4, fontsize=8)
plt.tight_layout()
plt.show()

imp = permutation_importance(
    gbr, test[cols], test["pv_actual"], n_repeats=5, random_state=SEED,
    scoring="neg_mean_absolute_error", n_jobs=4)
ranked = pd.Series(imp.importances_mean, index=cols).sort_values()

fig, ax = plt.subplots(figsize=(8, 6))
ranked.tail(15).plot.barh(ax=ax, color="tab:blue",
                          xerr=pd.Series(imp.importances_std, index=cols).reindex(
                              ranked.tail(15).index))
ax.set(title="Permutation importance, top 15", xlabel="increase in MAE [MW]",
       ylabel="")
plt.tight_layout()
plt.show()

WEATHER = ["shortwave_radiation", "direct_radiation", "diffuse_radiation",
           "temperature_2m", "cloud_cover", "cloud_idx"]

rng = np.random.default_rng(SEED)
day_mask = test["is_day"].to_numpy()
y_true = test["pv_actual"].to_numpy()
base = np.abs(gbr.predict(test[cols]).clip(0) - y_true)[day_mask].mean()

rows = []
for city in SITES:
    city_cols = [f"{w}_{city}" for w in WEATHER]
    ups = []
    for _ in range(5):
        X = test[cols].copy()
        X[city_cols] = X[city_cols].to_numpy()[rng.permutation(len(X))]
        ups.append(np.abs(gbr.predict(X).clip(0) - y_true)[day_mask].mean() - base)
    rows.append({"city": city, "dMAE_day": np.mean(ups), "std": np.std(ups)})

group = pd.DataFrame(rows).set_index("city").sort_values("dMAE_day", ascending=False)
print(f"daylight MAE {base:.0f} MW; increase when a city's weather is shuffled\n")
print(group.round(1))