#!.venv/bin/python """Backtest the local PV forecast model against measured PV production. Answers the question "is my PV forecast actually any good, and does a different configuration help?" without waiting for new forecasts to come true. Open-Meteo serves past weather in the same request as the forecast, so every variant can be scored right now against the meter readings EOS already holds. The comparison runs on past intervals, where the Open-Meteo rows are analysed rather than forecast weather. That isolates the error of the *PV model* (wrong peakpower, soiling, shading the horizon profile misses) from the error of the *weather forecast*. Only the former is systematic enough to fix by configuration. Requires: - ``general.latitude`` / ``general.longitude`` and ``pvforecast.planes`` configured - ``measurement.pv_production_emr_keys`` configured and fed with cumulative PV production meter readings [kWh] Usage: python scripts/pvforecast_backtest.py --days 30 python scripts/pvforecast_backtest.py --days 30 --tilt 88 --azimuth 175 """ import argparse import sys from pathlib import Path from typing import Any, Optional import numpy as np import pandas as pd # Add the src directory to sys.path so import akkudoktoreos works in all cases PROJECT_ROOT = Path(__file__).parent.parent SRC_DIR = PROJECT_ROOT / "src" sys.path.insert(0, str(SRC_DIR)) from akkudoktoreos.core.coreabc import get_config, get_measurement, singletons_init from akkudoktoreos.prediction.pvforecastakkudoktorlocal import ( PVForecastAkkudoktorLocal, PVForecastAkkudoktorLocalCommonSettings, ) from akkudoktoreos.utils.datetimeutil import to_datetime, to_duration ENSEMBLE = ["icon_seamless", "ecmwf_ifs025", "gfs_seamless"] # Variants scored against the meter. Each entry is a label plus the settings overrides # applied on top of the configured provider settings. VARIANTS: list[tuple[str, dict[str, Any]]] = [ ("best_match", {"weather_models": ["best_match"]}), ("ensemble", {"weather_models": ENSEMBLE}), ("ensemble + calibration", {"weather_models": ENSEMBLE, "calibration_enabled": True}), ( "ensemble + calibration (global only)", { "weather_models": ENSEMBLE, "calibration_enabled": True, "calibration_azimuth_bin_degrees": 0, }, ), ("ensemble, isotropic sky", {"weather_models": ENSEMBLE, "transposition_model": "isotropic"}), ("ensemble, no IAM", {"weather_models": ENSEMBLE, "apply_iam": False}), ] def score(modelled_kwh: np.ndarray, measured_kwh: np.ndarray) -> dict[str, float]: """Mean absolute error, bias and correlation over the daylight hours.""" error = modelled_kwh - measured_kwh measured_total = measured_kwh.sum() correlation = 0.0 if modelled_kwh.std() > 0 and measured_kwh.std() > 0: correlation = float(np.corrcoef(modelled_kwh, measured_kwh)[0, 1]) return { "mae": float(np.abs(error).mean()), "rmse": float(np.sqrt((error**2).mean())), "bias": float(error.mean()), "bias_pct": float(100.0 * error.sum() / measured_total) if measured_total > 0 else 0.0, "r": correlation, "model_kwh": float(modelled_kwh.sum()), "measured_kwh": float(measured_total), } def hourly_model(provider: PVForecastAkkudoktorLocal, data: Any, index: pd.DatetimeIndex) -> np.ndarray: """Run the chain and resample the AC power onto the measurement's hourly grid.""" frame = provider._forecast_frame(data) if frame.empty: return np.full(len(index), np.nan) # `ac_power` is a mean power per interval, so an hourly mean in W is Wh per hour. hourly = frame["ac_power"].resample("1h").mean().reindex(index) return hourly.to_numpy(dtype=float) / 1000.0 def main(days: int, tilt: Optional[float], azimuth: Optional[float]) -> int: singletons_init() config = get_config() measurement = get_measurement() if not config.measurement.pv_production_emr_keys: print( "measurement.pv_production_emr_keys is not configured - nothing to compare " "against. Configure it and feed cumulative PV production readings [kWh]." ) return 1 if not config.pvforecast.planes: print("pvforecast.planes is not configured.") return 1 if measurement.max_datetime is None: print("No measurements stored yet.") return 1 if tilt is not None: config.pvforecast.planes[0].surface_tilt = tilt if azimuth is not None: config.pvforecast.planes[0].surface_azimuth = azimuth end = measurement.max_datetime.start_of("hour") start = end.subtract(days=days) if start < measurement.min_datetime: start = measurement.min_datetime.start_of("hour").add(hours=1) measured_kwh = np.asarray( measurement.pv_production_total_kwh( start_datetime=start, end_datetime=end, interval=to_duration("1 hour") ), dtype=float, ) grid = pd.date_range( start=pd.Timestamp(start.in_timezone("UTC").isoformat()), periods=len(measured_kwh), freq="1h", ) print(f"Window: {start} .. {end} ({len(measured_kwh)} h)") print(f"Measured PV production: {measured_kwh.sum():.1f} kWh") provider = PVForecastAkkudoktorLocal(config=config, start_datetime=to_datetime()) baseline_settings = config.pvforecast.provider_settings.PVForecastAkkudoktorLocal if baseline_settings is None: baseline_settings = PVForecastAkkudoktorLocalCommonSettings() common = {"past_days": min(days + 1, 92), "calibration_days": days} rows = [] for label, overrides in VARIANTS: settings = baseline_settings.model_copy(update={**common, **overrides}) config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = settings try: data = provider._request_forecast(force_update=True) modelled_kwh = hourly_model(provider, data, grid) except Exception as exc: # noqa: BLE001 - one bad variant must not stop the rest print(f" {label}: failed ({exc})") continue # Only score hours where both sides exist and something was actually produced. usable = np.isfinite(modelled_kwh) & np.isfinite(measured_kwh) usable &= (modelled_kwh > 0.05) | (measured_kwh > 0.05) if usable.sum() < 12: print(f" {label}: too few usable hours ({int(usable.sum())})") continue rows.append((label, score(modelled_kwh[usable], measured_kwh[usable]), int(usable.sum()))) if not rows: print("No variant could be scored.") return 1 print(f"\n{'variant':38s} {'MAE':>7s} {'RMSE':>7s} {'bias':>8s} {'bias%':>7s} {'r':>6s}") print("-" * 78) for label, result, _ in sorted(rows, key=lambda row: row[1]["mae"]): print( f"{label:38s} {result['mae']:7.3f} {result['rmse']:7.3f} " f"{result['bias']:+8.3f} {result['bias_pct']:+6.1f}% {result['r']:6.3f}" ) print("\nMAE/RMSE/bias in kWh per hour. Lower MAE is better; bias% is the total") print("over- (+) or under-estimate (-) relative to the measured energy.") print( "\nNote: the calibrated variants are fitted on the same window they are scored\n" "on, so their advantage here is optimistic. Re-run with a longer --days to see\n" "how much of it survives." ) best_label, best, hours = min(rows, key=lambda row: row[1]["mae"]) print( f"\nBest: {best_label} - {best['model_kwh']:.1f} kWh modelled vs. " f"{best['measured_kwh']:.1f} kWh measured over {hours} scored hours." ) return 0 if __name__ == "__main__": parser = argparse.ArgumentParser(description="Backtest the local PV forecast model.") parser.add_argument("--days", type=int, default=30, help="Length of the window (default 30).") parser.add_argument("--tilt", type=float, default=None, help="Override plane 0 surface_tilt.") parser.add_argument( "--azimuth", type=float, default=None, help="Override plane 0 surface_azimuth." ) args = parser.parse_args() sys.exit(main(args.days, args.tilt, args.azimuth))