mirror of
https://github.com/Akkudoktor-EOS/EOS.git
synced 2026-10-09 07:56:40 +00:00
feat(pvforecast): local pvlib provider with measurement calibration
Add PVForecastAkkudoktorLocal, which runs the modelling chain inside EOS on raw Open-Meteo irradiance instead of calling a forecast service: solar position, horizon shading, plane transposition, incidence-angle modifier, cell temperature, PVWatts DC and inverter AC. It needs no API key and serves up to 16 days at 15-minute resolution from a single hourly request, which is what keeps `optimization.tail_horizon_hours` fed - services wrapping Open-Meteo cut the horizon much shorter. Several Open-Meteo models can be listed in `weather_models` and are averaged per variable at no extra request cost. With `calibration_enabled` the provider fits itself against `measurement.pv_production_emr_keys` over the past `calibration_days`: a global scale factor plus optional per-solar-azimuth factors, each weighted by modelled energy, shrunk toward the global factor by `calibration_prior_kwh` and clamped to `[calibration_min_factor, calibration_max_factor]`. The comparison runs on past intervals, where Open-Meteo serves analysed rather than forecast weather, so it corrects the error of the PV model and not that of the weather forecast. Calibration is a scale factor on the output and never touches `userhorizon`, `surface_tilt`, `surface_azimuth` or `peakpower`. The docs say so, and say why a short window and a plant fault inside it are the two ways to end up with a misleading factor. Also add `Measurement.pv_production_total_kwh()` alongside the existing load total, and `scripts/pvforecast_backtest.py`, which scores configuration variants against the stored meter readings without waiting for new forecasts to come true.
This commit is contained in:
Andreas
@@ -0,0 +1,195 @@
|
||||
#!.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))
|
||||
Reference in New Issue
Block a user