From f9763351220aea2a1bfc8226a5280399964f7ac2 Mon Sep 17 00:00:00 2001 From: Andreas Date: Sun, 6 Sep 2026 18:21:20 +0200 Subject: [PATCH] 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. --- CHANGELOG.md | 20 + docs/_generated/configpvforecast.md | 76 +- .../nodered_pv_measurement_function.js | 68 ++ docs/akkudoktoreos/prediction.md | 134 +++ scripts/pvforecast_backtest.py | 195 +++++ src/akkudoktoreos/measurement/measurement.py | 84 +- src/akkudoktoreos/prediction/prediction.py | 6 + src/akkudoktoreos/prediction/pvforecast.py | 6 + .../prediction/pvforecastakkudoktorlocal.py | 810 ++++++++++++++++++ .../prediction/weatheropenmeteo.py | 9 +- tests/test_prediction.py | 11 +- tests/test_pvforecastakkudoktorlocal.py | 291 +++++++ tests/test_weatheropenmeteo.py | 16 +- .../testdata/weatherforecast_openmeteo_1.json | 93 +- .../testdata/weatherforecast_openmeteo_2.json | 74 +- 15 files changed, 1814 insertions(+), 79 deletions(-) create mode 100644 docs/akkudoktoreos/nodered_pv_measurement_function.js create mode 100644 scripts/pvforecast_backtest.py create mode 100644 src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py create mode 100644 tests/test_pvforecastakkudoktorlocal.py diff --git a/CHANGELOG.md b/CHANGELOG.md index cf359711..10bc6b03 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -98,6 +98,26 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). day-ahead market prices are retained at their native hourly or quarter-hourly resolution and missing slots at the end of the optimization horizon are extended with weekly or daily seasonal ETS forecasts. A median fallback is used when the available history is too short for ETS. +- Add the `PVForecastAkkudoktorLocal` PV forecast provider, which runs the whole modelling chain + inside EOS with `pvlib` 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, serves up to 16 days at 15-minute resolution + from one hourly request - enough to feed `optimization.tail_horizon_hours` - and exposes albedo, + inverter efficiency and the module temperature coefficient as real configuration. Several + Open-Meteo models can be listed in `weather_models` and are averaged per variable at no extra + request cost. +- The local provider can calibrate itself against measured PV production (`calibration_enabled`). + It compares its own model against `measurement.pv_production_emr_keys` over the past + `calibration_days` and fits 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. Setting + `calibration_azimuth_bin_degrees` to 0 fits the global factor alone, which is what a short + window supports. +- Add `scripts/pvforecast_backtest.py`, which scores PV forecast configuration variants against + the stored meter readings straight away instead of waiting for new forecasts to come true, and + `Measurement.pv_production_total_kwh()` alongside the existing load total. ### Changed - Replace the fixed DEAP variation loop with adaptive genetic evolution. Crossover offspring may diff --git a/docs/_generated/configpvforecast.md b/docs/_generated/configpvforecast.md index a5776e00..f6d1811c 100644 --- a/docs/_generated/configpvforecast.md +++ b/docs/_generated/configpvforecast.md @@ -34,7 +34,8 @@ "PVForecastVrm": null, "PVForecastPVNode": null, "PVForecastForecastSolar": null, - "PVForecastSolcast": null + "PVForecastSolcast": null, + "PVForecastAkkudoktorLocal": null }, "planes": [ { @@ -102,7 +103,8 @@ "PVForecastVrm": null, "PVForecastPVNode": null, "PVForecastForecastSolar": null, - "PVForecastSolcast": null + "PVForecastSolcast": null, + "PVForecastAkkudoktorLocal": null }, "planes": [ { @@ -157,7 +159,8 @@ "PVForecastPVNode", "PVForecastForecastSolar", "PVForecastSolcast", - "PVForecastImport" + "PVForecastImport", + "PVForecastAkkudoktorLocal" ], "planes_peakpower": [ 5.0, @@ -192,6 +195,69 @@ ``` +### Common settings for the local (pvlib) PV forecast provider + + +:::{table} pvforecast::provider_settings::PVForecastAkkudoktorLocal +:widths: 10 10 5 5 30 +:align: left + +| Name | Type | Read-Only | Default | Description | +| ---- | ---- | --------- | ------- | ----------- | +| albedo | `float` | `rw` | `0.25` | Ground albedo used for planes that do not set their own. | +| apply_iam | `bool` | `rw` | `True` | Apply the ASHRAE incidence-angle modifier to the beam component. | +| calibration_azimuth_bin_degrees | `int` | `rw` | `15` | Width of the solar-azimuth bins for the correction. 0 fits a single global factor only. | +| calibration_days | `int` | `rw` | `30` | Length of the measurement window used to fit the correction. | +| calibration_enabled | `bool` | `rw` | `False` | Correct systematic model error against measured PV production. Requires `measurement.pv_production_emr_keys` to be configured and fed. Fits a global scale factor plus per-solar-azimuth factors, which is what catches near-field shading the horizon profile misses. | +| calibration_max_factor | `float` | `rw` | `1.5` | Upper clamp on any fitted correction factor. | +| calibration_min_factor | `float` | `rw` | `0.5` | Lower clamp on any fitted correction factor. | +| calibration_prior_kwh | `float` | `rw` | `5.0` | Shrinkage strength: a bin needs this much modelled energy before its own factor outweighs the global one. Higher is more conservative. | +| forecast_days | `Optional[int]` | `rw` | `None` | Forecast horizon in days (1-16). Leave empty to derive it from `prediction.hours`, which is what keeps the optimizer's tail horizon fed. | +| inverter_efficiency | `float` | `rw` | `0.96` | Nominal inverter efficiency (PVWatts eta_inv_nom). | +| past_days | `Optional[int]` | `rw` | `None` | Days of past data to request (0-92). Leave empty to derive it from `prediction.historic_hours`. | +| resolution_minutes | `int` | `rw` | `15` | Forecast resolution in minutes. 15 requests Open-Meteo's `minutely_15` block (natively resolved over Central Europe and North America, interpolated from hourly elsewhere); 60 requests the `hourly` block. | +| shift_to_interval_start | `bool` | `rw` | `True` | Open-Meteo stamps an interval mean with the interval END. EOS labels an interval by its START, so records are shifted back by one interval. Disable only to compare like-for-like against a provider that does not. | +| temperature_coefficient | `float` | `rw` | `-0.36` | Module power temperature coefficient in %/degC (negative). Matches the `cellCoEff` the akkudoktor.net forecast uses. | +| transposition_model | `str` | `rw` | `perez` | pvlib sky-diffuse transposition model: isotropic, klucher, haydavies, reindl, king or perez. | +| weather_models | `list[str]` | `rw` | `['best_match']` | Open-Meteo weather models to request. Listing more than one turns the input into a poor-man's ensemble: the members are averaged per variable, which is the cheapest reliable way to cut irradiance forecast error. Costs no extra API calls. | +::: + + + +**Example Input/Output** + + + +```json + { + "pvforecast": { + "provider_settings": { + "PVForecastAkkudoktorLocal": { + "resolution_minutes": 15, + "forecast_days": null, + "past_days": null, + "weather_models": [ + "best_match" + ], + "transposition_model": "perez", + "albedo": 0.25, + "inverter_efficiency": 0.96, + "temperature_coefficient": -0.36, + "apply_iam": true, + "shift_to_interval_start": true, + "calibration_enabled": true, + "calibration_days": 30, + "calibration_azimuth_bin_degrees": 15, + "calibration_prior_kwh": 5.0, + "calibration_min_factor": 0.5, + "calibration_max_factor": 1.5 + } + } + } + } +``` + + ### Common settings for the Solcast PV forecast provider @@ -368,6 +434,7 @@ | ---- | ---- | --------- | ------- | ----------- | | PVForecastForecastSolar | `Optional[akkudoktoreos.prediction.pvforecastforecastsolar.PVForecastForecastSolarCommonSettings]` | `rw` | `None` | PVForecastForecastSolar settings | | PVForecastImport | `Optional[akkudoktoreos.prediction.pvforecastimport.PVForecastImportCommonSettings]` | `rw` | `None` | PVForecastImport settings | +| PVForecastAkkudoktorLocal | `Optional[akkudoktoreos.prediction.pvforecastlocal.PVForecastAkkudoktorLocalCommonSettings]` | `rw` | `None` | PVForecastAkkudoktorLocal settings | | PVForecastPVNode | `Optional[akkudoktoreos.prediction.pvforecastpvnode.PVForecastPVNodeCommonSettings]` | `rw` | `None` | PVForecastPVNode settings | | PVForecastSolcast | `Optional[akkudoktoreos.prediction.pvforecastsolcast.PVForecastSolcastCommonSettings]` | `rw` | `None` | PVForecastSolcast settings | | PVForecastVrm | `Optional[akkudoktoreos.prediction.pvforecastvrm.PVForecastVrmCommonSettings]` | `rw` | `None` | PVForecastVrm settings | @@ -387,7 +454,8 @@ "PVForecastVrm": null, "PVForecastPVNode": null, "PVForecastForecastSolar": null, - "PVForecastSolcast": null + "PVForecastSolcast": null, + "PVForecastAkkudoktorLocal": null } } } diff --git a/docs/akkudoktoreos/nodered_pv_measurement_function.js b/docs/akkudoktoreos/nodered_pv_measurement_function.js new file mode 100644 index 00000000..80ee81d6 --- /dev/null +++ b/docs/akkudoktoreos/nodered_pv_measurement_function.js @@ -0,0 +1,68 @@ +// Node-RED function: cumulative PV production meter readings -> EOS measurement API. +// +// Strang: +// inject (repeat 3600 s, msg.topic = die SQL aus dem Tutorial) +// -> mysql "MariaDB" (DB `sensor`) +// -> DIESE function +// -> http request (Method: "- set by msg.method -") +// +// EOS speichert unter `pv_production_emr_keys` Zaehlerstaende (EMR) in kWh und +// bildet die Differenzen selbst. Aus den Momentanleistungen in `data.solarallpower` +// muss also erst ein monoton steigender Zaehler werden - das macht die SQL. +// +// Wichtig: das Startdatum in der SQL bleibt FEST. Ein rollendes +// `NOW() - INTERVAL n DAY` verschiebt den Nullpunkt der kumulativen Summe bei +// jedem Lauf, und EOS liest den Sprung als Produktion. + +const BASE_URL = "http://192.168.1.151:8503"; +const KEY = "pv_produktion_emr"; +const TZ = "Europe/Berlin"; + +function toIsoWithOffset(value) { + // DATE_FORMAT() kommt als String in lokaler Zeit zurueck. Den Offset aus dem + // Datum selbst bilden, damit CEST und CET beide stimmen. + const date = value instanceof Date ? value : new Date(String(value).replace(" ", "T")); + const pad = n => String(Math.trunc(Math.abs(n))).padStart(2, "0"); + const offsetMin = -date.getTimezoneOffset(); + const sign = offsetMin >= 0 ? "+" : "-"; + return `${date.getFullYear()}-${pad(date.getMonth() + 1)}-${pad(date.getDate())}` + + `T${pad(date.getHours())}:${pad(date.getMinutes())}:${pad(date.getSeconds())}` + + `${sign}${pad(offsetMin / 60)}:${pad(offsetMin % 60)}`; +} + +const rows = Array.isArray(msg.payload) ? msg.payload : []; +const data = {}; +let last = null; +let skipped = 0; +for (const row of rows) { + const value = Number(row.emr); + if (!Number.isFinite(value)) { + skipped += 1; + continue; + } + // Ein Zaehler laeuft nur vorwaerts. Ein Rueckschritt bedeutet eine Luecke oder + // einen kaputten Messwert - EOS wuerde daraus eine negative Produktionsstunde + // machen. + if (last !== null && value < last) { + skipped += 1; + continue; + } + last = value; + data[toIsoWithOffset(row.ts)] = Number(value.toFixed(6)); +} + +const count = Object.keys(data).length; +if (count === 0) { + node.warn("Keine PV-Messwerte gefunden - topic und Zeitfenster in der SQL pruefen."); + return null; +} +if (skipped > 0) { + node.warn(`${skipped} Zeilen uebersprungen (nicht-numerisch oder Zaehler rueckwaerts).`); +} +node.status({ text: `${count} EMR-Werte, letzter ${last.toFixed(1)} kWh` }); + +msg.method = "PUT"; +msg.url = `${BASE_URL}/v1/measurement/series?key=${encodeURIComponent(KEY)}`; +msg.headers = { "Content-Type": "application/json" }; +msg.payload = { data: data, dtype: "float64", tz: TZ }; +return msg; diff --git a/docs/akkudoktoreos/prediction.md b/docs/akkudoktoreos/prediction.md index 3d15269f..7c0ef7e3 100644 --- a/docs/akkudoktoreos/prediction.md +++ b/docs/akkudoktoreos/prediction.md @@ -489,6 +489,7 @@ Configuration options: - `PVForecastForecastSolar`: Retrieves forecasts from the free Forecast.Solar API. - `PVForecastSolcast`: Retrieves forecasts from the Solcast rooftop-site API. - `PVForecastImport`: Imports from a file or JSON string or by endpoint data provision. + - `PVForecastAkkudoktorLocal`: Computes the forecast inside EOS from Open-Meteo weather with pvlib. - `planes[].surface_tilt`: Tilt angle from horizontal plane. Ignored for two-axis tracking. - `planes[].surface_azimuth`: Orientation (azimuth angle) of the (fixed) plane. @@ -522,6 +523,17 @@ Configuration options: - `provider_settings.PVForecastForecastSolar.api_key`: Forecast.Solar API key (optional). - `provider_settings.PVForecastSolcast.api_key`: Solcast API key. - `provider_settings.PVForecastSolcast.site_id`: Solcast rooftop resource (site) id. + - `provider_settings.PVForecastAkkudoktorLocal.resolution_minutes`: 15 or 60. + - `provider_settings.PVForecastAkkudoktorLocal.forecast_days`: 1-16, empty derives it from `prediction.hours`. + - `provider_settings.PVForecastAkkudoktorLocal.past_days`: 0-92, empty derives it from `prediction.historic_hours`. + - `provider_settings.PVForecastAkkudoktorLocal.weather_models`: Open-Meteo models; several are averaged. + - `provider_settings.PVForecastAkkudoktorLocal.transposition_model`: pvlib sky-diffuse model. + - `provider_settings.PVForecastAkkudoktorLocal.albedo`: Fallback albedo for planes without their own. + - `provider_settings.PVForecastAkkudoktorLocal.inverter_efficiency`: Nominal inverter efficiency. + - `provider_settings.PVForecastAkkudoktorLocal.temperature_coefficient`: Module power coefficient in %/degC. + - `provider_settings.PVForecastAkkudoktorLocal.apply_iam`: Apply the ASHRAE incidence-angle modifier. + - `provider_settings.PVForecastAkkudoktorLocal.shift_to_interval_start`: Relabel Open-Meteo interval-end stamps. + - `provider_settings.PVForecastAkkudoktorLocal.calibration_*`: Self-calibration against measured PV production. --- @@ -752,6 +764,128 @@ The prediction keys for the PV forecast data are: - `pvforecast_ac_power`: Total AC power (W). - `pvforecast_dc_power`: Total DC power (W). +### PVForecastAkkudoktorLocal Provider + +The `PVForecastAkkudoktorLocal` provider does not call a PV forecast service at all. It fetches raw +irradiance and weather from [Open-Meteo](https://open-meteo.com) and runs the whole modelling +chain locally with `pvlib`: + + solar position -> horizon shading -> transposition to the module plane -> + incidence-angle modifier -> cell temperature -> PVWatts DC -> inverter AC + +Three properties make it the right default for long-horizon optimization: + +- **Horizon.** Up to 16 forecast days at 15-minute resolution from a single request. Services + that wrap Open-Meteo cut the horizon much shorter, which starves + `optimization.tail_horizon_hours`. +- **Call budget.** One request per hour against a ~10k/day non-commercial budget, instead of + competing for someone else's upstream quota. +- **Honest parameters.** `albedo`, inverter efficiency and the module temperature coefficient are + real configuration rather than constants baked into a service URL. + +No API key is required. The location comes from `general.latitude`/`longitude` and the geometry +from the configured `planes`, including `userhorizon`, `trackingtype`, `mountingplace` and `loss`. + +```python + { + "pvforecast": { + "provider": "PVForecastAkkudoktorLocal", + "provider_settings": { + "PVForecastAkkudoktorLocal": { + "resolution_minutes": 15, + "weather_models": ["icon_seamless", "ecmwf_ifs025", "gfs_seamless"], + "calibration_enabled": true + } + } + } + } +``` + +#### Improving accuracy + +**Model ensemble.** Listing several models in `weather_models` averages them per variable. This +is the cheapest reliable way to cut irradiance forecast error and costs no extra API calls, +because Open-Meteo returns all members in the same response. `icon_seamless` (DWD, strong over +Central Europe), `ecmwf_ifs025` and `gfs_seamless` are a reasonable trio. + +**Self-calibration.** With `calibration_enabled` the provider compares its own model against +measured PV production over the past `calibration_days` and fits a correction: + +- a **global scale factor**, which absorbs a wrong `peakpower`, soiling, degradation and any + systematic offset in the loss assumption; +- **per-solar-azimuth factors**, which absorb near-field shading that a coarse `userhorizon` + cannot express - a chimney, a tree, a neighbouring roof. + +Each bin is weighted by its modelled energy and shrunk toward the global factor by +`calibration_prior_kwh`, so a thinly sampled bin cannot swing the forecast on its own, and every +factor is clamped to `[calibration_min_factor, calibration_max_factor]` so a broken meter cannot +either. The fitted factors and the resulting change in mean absolute error are logged at INFO +level on every update. + +Calibration requires `measurement.pv_production_emr_keys` to be configured and fed with +cumulative PV production meter readings in kWh: + +```python + { + "measurement": { + "pv_production_emr_keys": ["pv1_emr"] + } + } +``` + +Note what this does and does not correct. The comparison runs on past intervals, where the +Open-Meteo rows are analysed rather than forecast weather, so it isolates the error of the *PV +model* from the error of the *weather forecast*. That is deliberate: only the former is +systematic enough to correct. A cloudy day that the weather model got wrong stays wrong. + +**Choosing the window.** `calibration_days` trades responsiveness against stability. A long +window averages more weather and gives a steadier factor, but it also reaches back into the +plant's own history - and calibration cannot tell a modelling error from a *plant* error. A +window that spans a string outage, an inverter derating or a period of heavy soiling will fit +that fault as if it were a permanent property of the installation, and the forecast stays +depressed long after the plant recovered. Before trusting a factor, compare modelled against +measured energy *per day*: a run of days at a markedly different ratio is a plant event, and +the window should start after it. + +**Bins need data.** The per-azimuth factors are only worth fitting when the window holds enough +daylight hours to populate the bins - roughly a few hundred, so several weeks at the default +15 degrees. With a short window, set `calibration_azimuth_bin_degrees` to 0 to fit the global +factor alone. One well-determined number beats twenty-four noisy ones. + +**Geometry is out of scope.** Calibration is a scale factor on the model's output; it never +touches `userhorizon`, `surface_tilt`, `surface_azimuth` or `peakpower`. That makes it the right +tool for *multiplicative* errors and the wrong one for geometric errors. Horizon shading in +particular gates the beam component as a hard function of both solar azimuth *and* elevation, so +a per-azimuth scale factor cannot move the edge of the shadow to where it belongs, and what it +learns in one season is wrong in the next, when the sun crosses the same azimuth at a different +height. A wrong horizon should be corrected in `userhorizon`, not calibrated away. + +`scripts/pvforecast_backtest.py` scores configuration variants against the stored meter readings +without waiting for new forecasts to come true, which is the quickest way to test a geometry +change: + +```bash + python scripts/pvforecast_backtest.py --days 30 --tilt 88 --azimuth 175 +``` + +#### Conventions + +Two timing conventions are handled explicitly and are worth knowing when comparing against other +providers: + +- Open-Meteo radiation values are the mean over the **preceding** interval, so the representative + sun position for a value stamped `t` is `t - interval/2`. +- EOS records label an interval by its **start**, so a value stamped `t` by Open-Meteo is stored + at `t - interval`. Set `shift_to_interval_start` to false to keep the raw stamps. + +Note also that Open-Meteo's `direct_radiation` is beam irradiance on the *horizontal* plane; the +DNI this chain needs is `direct_normal_irradiance`. + +The prediction keys for the PV forecast data are: + +- `pvforecast_ac_power`: Total AC power (W). +- `pvforecast_dc_power`: Total DC power (W). + ### PVForecastForecastSolar Provider The `PVForecastForecastSolar` provider retrieves PV power forecasts from the free diff --git a/scripts/pvforecast_backtest.py b/scripts/pvforecast_backtest.py new file mode 100644 index 00000000..02f06aa0 --- /dev/null +++ b/scripts/pvforecast_backtest.py @@ -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)) diff --git a/src/akkudoktoreos/measurement/measurement.py b/src/akkudoktoreos/measurement/measurement.py index dbbdca41..3c6d1328 100644 --- a/src/akkudoktoreos/measurement/measurement.py +++ b/src/akkudoktoreos/measurement/measurement.py @@ -203,22 +203,25 @@ class Measurement(SingletonMixin, DataImportMixin, DataSequence): logger.debug(debug_msg) return energy_array - def load_total_kwh( + def _total_kwh( self, + emr_keys: Optional[list[str]], + label: str, start_datetime: Optional[DateTime] = None, end_datetime: Optional[DateTime] = None, interval: Optional[Duration] = None, ) -> NDArray[Shape["*"], Any]: - """Calculate a total load energy values array indexed by fixed time intervals from load metering data within an optional date range. + """Sum the per-interval energy of several meter reading keys. Args: - start_datetime (datetime, optional): The start date for filtering the load data (inclusive). - end_datetime (datetime, optional): The end date for filtering the load data (exclusive). + emr_keys: The configured energy meter reading keys to sum, or None. + label: Name of the summed quantity, used for debug logging only. + start_datetime (datetime, optional): The start date for filtering (inclusive). + end_datetime (datetime, optional): The end date for filtering (exclusive). interval (duration, optional): The fixed time interval. Defaults to 1 hour. Returns: - np.ndarray: A NumPy Array of the total load energy [kWh] per interval values calculated from - the load meter readings. + np.ndarray: A NumPy Array of the total energy [kWh] per interval. """ if interval is None: interval = to_duration("1 hour") @@ -236,24 +239,73 @@ class Measurement(SingletonMixin, DataImportMixin, DataSequence): if end_datetime is None: end_datetime = self.max_datetime.add(seconds=1) size = self._interval_count(start_datetime, end_datetime, interval) - load_total_kwh_array = np.zeros(size) + total_kwh_array = np.zeros(size) - # Loop through all loads - if isinstance(self.config.measurement.load_emr_keys, list): - for key in self.config.measurement.load_emr_keys: - # Calculate load per interval - load_array = self._energy_from_meter_readings( + if isinstance(emr_keys, list): + for key in emr_keys: + # Calculate energy per interval + energy_array = self._energy_from_meter_readings( key=key, start_datetime=start_datetime, end_datetime=end_datetime, interval=interval, ) - # Add calculated load to total load - load_total_kwh_array += load_array - debug_msg = f"Total load '{key}' calculation: {load_total_kwh_array}" + # Add to the total + total_kwh_array += energy_array + debug_msg = f"Total {label} '{key}' calculation: {total_kwh_array}" logger.debug(debug_msg) - return load_total_kwh_array + return total_kwh_array + + def load_total_kwh( + self, + start_datetime: Optional[DateTime] = None, + end_datetime: Optional[DateTime] = None, + interval: Optional[Duration] = None, + ) -> NDArray[Shape["*"], Any]: + """Calculate a total load energy values array indexed by fixed time intervals from load metering data within an optional date range. + + Args: + start_datetime (datetime, optional): The start date for filtering the load data (inclusive). + end_datetime (datetime, optional): The end date for filtering the load data (exclusive). + interval (duration, optional): The fixed time interval. Defaults to 1 hour. + + Returns: + np.ndarray: A NumPy Array of the total load energy [kWh] per interval values calculated from + the load meter readings. + """ + return self._total_kwh( + emr_keys=self.config.measurement.load_emr_keys, + label="load", + start_datetime=start_datetime, + end_datetime=end_datetime, + interval=interval, + ) + + def pv_production_total_kwh( + self, + start_datetime: Optional[DateTime] = None, + end_datetime: Optional[DateTime] = None, + interval: Optional[Duration] = None, + ) -> NDArray[Shape["*"], Any]: + """Calculate total PV production energy per interval from PV meter readings. + + Args: + start_datetime (datetime, optional): The start date for filtering the data (inclusive). + end_datetime (datetime, optional): The end date for filtering the data (exclusive). + interval (duration, optional): The fixed time interval. Defaults to 1 hour. + + Returns: + np.ndarray: A NumPy Array of the total PV production energy [kWh] per interval + calculated from the PV production meter readings. + """ + return self._total_kwh( + emr_keys=self.config.measurement.pv_production_emr_keys, + label="PV production", + start_datetime=start_datetime, + end_datetime=end_datetime, + interval=interval, + ) # ----------------------- Measurement Database Protocol --------------------- diff --git a/src/akkudoktoreos/prediction/prediction.py b/src/akkudoktoreos/prediction/prediction.py index 39073ee9..b35b6352 100644 --- a/src/akkudoktoreos/prediction/prediction.py +++ b/src/akkudoktoreos/prediction/prediction.py @@ -53,6 +53,7 @@ from akkudoktoreos.prediction.predictionabc import PredictionContainer from akkudoktoreos.prediction.pvforecastakkudoktor import PVForecastAkkudoktor from akkudoktoreos.prediction.pvforecastforecastsolar import PVForecastForecastSolar from akkudoktoreos.prediction.pvforecastimport import PVForecastImport +from akkudoktoreos.prediction.pvforecastakkudoktorlocal import PVForecastAkkudoktorLocal from akkudoktoreos.prediction.pvforecastpvnode import PVForecastPVNode from akkudoktoreos.prediction.pvforecastsolcast import PVForecastSolcast from akkudoktoreos.prediction.pvforecastvrm import PVForecastVrm @@ -103,6 +104,7 @@ pvforecast_pvnode = PVForecastPVNode() pvforecast_forecastsolar = PVForecastForecastSolar() pvforecast_solcast = PVForecastSolcast() pvforecast_import = PVForecastImport() +pvforecast_akkudoktor_local = PVForecastAkkudoktorLocal() weather_brightsky = WeatherBrightSky() weather_clearoutside = WeatherClearOutside() weather_openmeteo = WeatherOpenMeteo() @@ -134,6 +136,7 @@ def prediction_providers() -> ( PVForecastForecastSolar, PVForecastSolcast, PVForecastImport, + PVForecastAkkudoktorLocal, WeatherBrightSky, WeatherClearOutside, WeatherOpenMeteo, @@ -168,6 +171,7 @@ def prediction_providers() -> ( pvforecast_forecastsolar, \ pvforecast_solcast, \ pvforecast_import, \ + pvforecast_akkudoktor_local, \ weather_brightsky, \ weather_clearoutside, \ weather_openmeteo, \ @@ -197,6 +201,7 @@ def prediction_providers() -> ( pvforecast_forecastsolar, pvforecast_solcast, pvforecast_import, + pvforecast_akkudoktor_local, weather_brightsky, weather_clearoutside, weather_openmeteo, @@ -231,6 +236,7 @@ class Prediction(PredictionContainer): PVForecastForecastSolar, PVForecastSolcast, PVForecastImport, + PVForecastAkkudoktorLocal, WeatherBrightSky, WeatherClearOutside, WeatherOpenMeteo, diff --git a/src/akkudoktoreos/prediction/pvforecast.py b/src/akkudoktoreos/prediction/pvforecast.py index 61381e73..292cb0e8 100644 --- a/src/akkudoktoreos/prediction/pvforecast.py +++ b/src/akkudoktoreos/prediction/pvforecast.py @@ -11,6 +11,7 @@ from akkudoktoreos.prediction.pvforecastforecastsolar import ( PVForecastForecastSolarCommonSettings, ) from akkudoktoreos.prediction.pvforecastimport import PVForecastImportCommonSettings +from akkudoktoreos.prediction.pvforecastakkudoktorlocal import PVForecastAkkudoktorLocalCommonSettings from akkudoktoreos.prediction.pvforecastpvnode import PVForecastPVNodeCommonSettings from akkudoktoreos.prediction.pvforecastsolcast import PVForecastSolcastCommonSettings from akkudoktoreos.prediction.pvforecastvrm import PVForecastVrmCommonSettings @@ -30,6 +31,7 @@ def pvforecast_provider_ids() -> list[str]: "PVForecastPVNode", "PVForecastForecastSolar", "PVForecastSolcast", + "PVForecastAkkudoktorLocal", ] return [ @@ -203,6 +205,10 @@ class PVForecastCommonProviderSettings(SettingsBaseModel): default=None, json_schema_extra={"description": "PVForecastSolcast settings", "examples": [None]}, ) + PVForecastAkkudoktorLocal: Optional[PVForecastAkkudoktorLocalCommonSettings] = Field( + default=None, + json_schema_extra={"description": "PVForecastAkkudoktorLocal settings", "examples": [None]}, + ) class PVForecastCommonSettings(SettingsBaseModel): diff --git a/src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py b/src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py new file mode 100644 index 00000000..d31b2e55 --- /dev/null +++ b/src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py @@ -0,0 +1,810 @@ +"""Native PV power forecast computed inside EOS from Open-Meteo weather with pvlib. + +Unlike the other PV forecast providers this one does not ask a third-party service +for PV power. It fetches raw irradiance and weather from Open-Meteo and runs the +whole modelling chain locally with pvlib: + + solar position -> horizon shading -> transposition to the module plane -> + incidence-angle modifier -> cell temperature -> PVWatts DC -> inverter AC + +Why this exists: + +* **Horizon.** Open-Meteo serves up to 16 forecast days at 15-minute resolution in a + single request. Services that wrap it (including api.akkudoktor.net) cut the + horizon much shorter, which starves ``optimization.tail_horizon_hours``. +* **Call budget.** One request per hour against a ~10k/day non-commercial budget, + instead of competing for someone else's upstream quota. +* **Honest parameters.** ``albedo``, inverter efficiency and the module temperature + coefficient become real configuration instead of constants baked into a URL. + +Two conventions matter and are handled explicitly here: + +* Open-Meteo radiation values are the mean over the **preceding** interval, so the + representative sun position for a value stamped ``t`` is ``t - interval/2``. +* EOS records label an interval by its **start** (``key_to_array`` resamples with + left-labelled buckets), so a value stamped ``t`` by Open-Meteo is stored at + ``t - interval``. Set ``shift_to_interval_start`` to False to keep the raw stamps. + +Note also that ``direct_radiation`` in the Open-Meteo API is beam irradiance on the +*horizontal* plane; the DNI this chain needs is ``direct_normal_irradiance``. +""" + +import math +from typing import Any, Optional + +import numpy as np +import pandas as pd +import pendulum +import pvlib +import requests +from loguru import logger +from pydantic import Field, field_validator + +from akkudoktoreos.config.configabc import SettingsBaseModel +from akkudoktoreos.core.cache import cache_in_file +from akkudoktoreos.prediction.pvforecastabc import PVForecastProvider +from akkudoktoreos.utils.datetimeutil import compare_datetimes, to_datetime, to_duration + +OPENMETEO_URL = "https://api.open-meteo.com/v1/forecast" + +# Open-Meteo variables the pvlib chain needs. `direct_normal_irradiance` is the DNI; +# `direct_radiation` (beam on the horizontal) would be wrong here. +OPENMETEO_VARIABLES = ( + "temperature_2m", + "relative_humidity_2m", + "wind_speed_10m", + "shortwave_radiation", + "diffuse_radiation", + "direct_normal_irradiance", +) + +# Open-Meteo hard limits for a single forecast request. +MAX_FORECAST_DAYS = 16 +MAX_PAST_DAYS = 92 + +TRANSPOSITION_MODELS = ("isotropic", "klucher", "haydavies", "reindl", "king", "perez") + +# pvlib SAPM cell temperature parameter set per EOS `mountingplace`. +MOUNTING_TEMPERATURE_MODEL = { + "free": "open_rack_glass_glass", + "building": "close_mount_glass_glass", +} + + +class PVForecastAkkudoktorLocalCommonSettings(SettingsBaseModel): + """Common settings for the local (pvlib) PV forecast provider.""" + + resolution_minutes: int = Field( + default=15, + json_schema_extra={ + "description": ( + "Forecast resolution in minutes. 15 requests Open-Meteo's `minutely_15` " + "block (natively resolved over Central Europe and North America, " + "interpolated from hourly elsewhere); 60 requests the `hourly` block." + ), + "examples": [15, 60], + }, + ) + forecast_days: Optional[int] = Field( + default=None, + ge=1, + le=MAX_FORECAST_DAYS, + json_schema_extra={ + "description": ( + "Forecast horizon in days (1-16). Leave empty to derive it from " + "`prediction.hours`, which is what keeps the optimizer's tail horizon fed." + ), + "examples": [None, 7], + }, + ) + past_days: Optional[int] = Field( + default=None, + ge=0, + le=MAX_PAST_DAYS, + json_schema_extra={ + "description": ( + "Days of past data to request (0-92). Leave empty to derive it from " + "`prediction.historic_hours`." + ), + "examples": [None, 3], + }, + ) + weather_models: list[str] = Field( + default=["best_match"], + min_length=1, + json_schema_extra={ + "description": ( + "Open-Meteo weather models to request. Listing more than one turns the " + "input into a poor-man's ensemble: the members are averaged per variable, " + "which is the cheapest reliable way to cut irradiance forecast error. " + "Costs no extra API calls." + ), + "examples": [ + ["best_match"], + ["icon_seamless", "ecmwf_ifs025", "gfs_seamless"], + ], + }, + ) + transposition_model: str = Field( + default="perez", + json_schema_extra={ + "description": ( + "pvlib sky-diffuse transposition model: isotropic, klucher, haydavies, " + "reindl, king or perez." + ), + "examples": ["perez", "haydavies"], + }, + ) + albedo: float = Field( + default=0.25, + ge=0.0, + le=1.0, + json_schema_extra={ + "description": "Ground albedo used for planes that do not set their own.", + "examples": [0.25, 0.2], + }, + ) + inverter_efficiency: float = Field( + default=0.96, + gt=0.0, + le=1.0, + json_schema_extra={ + "description": "Nominal inverter efficiency (PVWatts eta_inv_nom).", + "examples": [0.96, 0.94], + }, + ) + temperature_coefficient: float = Field( + default=-0.36, + json_schema_extra={ + "description": ( + "Module power temperature coefficient in %/degC (negative). Matches the " + "`cellCoEff` the akkudoktor.net forecast uses." + ), + "examples": [-0.36, -0.29], + }, + ) + apply_iam: bool = Field( + default=True, + json_schema_extra={ + "description": "Apply the ASHRAE incidence-angle modifier to the beam component.", + "examples": [True], + }, + ) + shift_to_interval_start: bool = Field( + default=True, + json_schema_extra={ + "description": ( + "Open-Meteo stamps an interval mean with the interval END. EOS labels an " + "interval by its START, so records are shifted back by one interval. " + "Disable only to compare like-for-like against a provider that does not." + ), + "examples": [True], + }, + ) + + calibration_enabled: bool = Field( + default=False, + json_schema_extra={ + "description": ( + "Correct systematic model error against measured PV production. Requires " + "`measurement.pv_production_emr_keys` to be configured and fed. Fits a " + "global scale factor plus per-solar-azimuth factors, which is what catches " + "near-field shading the horizon profile misses." + ), + "examples": [True], + }, + ) + calibration_days: int = Field( + default=30, + ge=1, + le=MAX_PAST_DAYS, + json_schema_extra={ + "description": "Length of the measurement window used to fit the correction.", + "examples": [30, 14], + }, + ) + calibration_azimuth_bin_degrees: int = Field( + default=15, + ge=0, + le=180, + json_schema_extra={ + "description": ( + "Width of the solar-azimuth bins for the correction. 0 fits a single " + "global factor only." + ), + "examples": [15, 30, 0], + }, + ) + calibration_prior_kwh: float = Field( + default=5.0, + ge=0.0, + json_schema_extra={ + "description": ( + "Shrinkage strength: a bin needs this much modelled energy before its own " + "factor outweighs the global one. Higher is more conservative." + ), + "examples": [5.0, 20.0], + }, + ) + calibration_min_factor: float = Field( + default=0.5, + gt=0.0, + json_schema_extra={ + "description": "Lower clamp on any fitted correction factor.", + "examples": [0.5], + }, + ) + calibration_max_factor: float = Field( + default=1.5, + gt=0.0, + json_schema_extra={ + "description": "Upper clamp on any fitted correction factor.", + "examples": [1.5], + }, + ) + + @field_validator("resolution_minutes") + @classmethod + def validate_resolution(cls, value: int) -> int: + if value not in (15, 60): + raise ValueError(f"resolution_minutes must be 15 or 60, got {value}") + return value + + @field_validator("transposition_model") + @classmethod + def validate_transposition_model(cls, value: str) -> str: + if value not in TRANSPOSITION_MODELS: + raise ValueError( + f"Invalid transposition_model '{value}', expected one of {TRANSPOSITION_MODELS}" + ) + return value + + +class PVForecastAkkudoktorLocal(PVForecastProvider): + """Compute the PV forecast locally from Open-Meteo irradiance using pvlib.""" + + @classmethod + def provider_id(cls) -> str: + """Return the unique identifier for the PV-Forecast-Provider.""" + return "PVForecastAkkudoktorLocal" + + @property + def _settings(self) -> PVForecastAkkudoktorLocalCommonSettings: + settings = self.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal + if settings is None: + settings = PVForecastAkkudoktorLocalCommonSettings() + return settings + + # ------------------------------------------------------------------ request + + def _horizon_days(self) -> tuple[int, int]: + """Resolve (forecast_days, past_days), deriving them from prediction config.""" + settings = self._settings + + forecast_days = settings.forecast_days + if forecast_days is None: + # +1 day so the last requested hour is still covered when the run starts + # late in the day. + hours = self.config.prediction.hours or 48 + forecast_days = min(MAX_FORECAST_DAYS, max(1, math.ceil(hours / 24) + 1)) + + past_days = settings.past_days + if past_days is None: + historic_hours = self.config.prediction.historic_hours or 0 + past_days = math.ceil(historic_hours / 24) + if settings.calibration_enabled: + # The fit compares modelled against measured power over the same past + # intervals, so the weather for that window has to come back with the request. + past_days = max(past_days, settings.calibration_days) + + return forecast_days, min(MAX_PAST_DAYS, past_days) + + @cache_in_file(with_ttl="1 hour") + def _request_forecast(self) -> Any: + """Fetch raw irradiance and weather from Open-Meteo.""" + latitude = self.config.general.latitude + longitude = self.config.general.longitude + if latitude is None or longitude is None: + raise ValueError("PVForecastAkkudoktorLocal needs general.latitude and general.longitude") + + settings = self._settings + block = "minutely_15" if settings.resolution_minutes == 15 else "hourly" + forecast_days, past_days = self._horizon_days() + + params = { + "latitude": latitude, + "longitude": longitude, + block: ",".join(OPENMETEO_VARIABLES), + # Ask for UTC so the returned stamps are unambiguous across DST changes. + "timezone": "UTC", + # pvlib's SAPM cell temperature model wants m/s; Open-Meteo defaults to km/h. + "wind_speed_unit": "ms", + "forecast_days": forecast_days, + "past_days": past_days, + "models": ",".join(settings.weather_models), + } + + try: + response = requests.get(OPENMETEO_URL, params=params, timeout=30) + logger.debug(f"Requesting Open-Meteo forecast: {response.url}") + response.raise_for_status() + except requests.RequestException as e: + logger.error(f"Failed to fetch weather for local pvforecast: {e}") + raise RuntimeError("Failed to fetch weather from Open-Meteo API") from e + + data = response.json() + if block not in data: + raise ValueError( + f"Open-Meteo response is missing the '{block}' block: {list(data.keys())}" + ) + + self.update_datetime = to_datetime(in_timezone=self.config.general.timezone) + return data + + # ------------------------------------------------------------------ weather + + def _weather_frame(self, data: Any) -> pd.DataFrame: + """Build a UTC-indexed weather frame from the Open-Meteo response.""" + block = "minutely_15" if self._settings.resolution_minutes == 15 else "hourly" + raw = data[block] + + index = pd.DatetimeIndex( + [pendulum.parse(str(t), tz="UTC") for t in raw["time"]], name="time" + ).tz_convert("UTC") + + frame = pd.DataFrame(index=index) + for variable in OPENMETEO_VARIABLES: + # With a single model Open-Meteo returns the bare variable name; with several + # it suffixes each member (`shortwave_radiation_icon_seamless`). Average the + # members that actually came back - not every model carries every variable. + members = [ + pd.to_numeric(pd.Series(raw[key], index=index), errors="coerce") + for key in raw + if key == variable or key.startswith(f"{variable}_") + ] + if not members: + raise ValueError(f"Open-Meteo response is missing '{variable}'") + frame[variable] = pd.concat(members, axis=1).mean(axis=1, skipna=True) + + # Night rows and occasional gaps arrive as null. Irradiance is genuinely zero + # then; temperature and wind are interpolated so the temperature model stays + # defined instead of poisoning the whole row with NaN. + for variable in ("shortwave_radiation", "diffuse_radiation", "direct_normal_irradiance"): + frame[variable] = frame[variable].fillna(0.0).clip(lower=0.0) + for variable in ("temperature_2m", "relative_humidity_2m", "wind_speed_10m"): + frame[variable] = frame[variable].interpolate(method="time").ffill().bfill() + frame["wind_speed_10m"] = frame["wind_speed_10m"].fillna(1.0).clip(lower=0.0) + frame["temperature_2m"] = frame["temperature_2m"].fillna(15.0) + + return frame + + # ------------------------------------------------------------------ geometry + + @staticmethod + def _horizon_elevation(userhorizon: list[float], solar_azimuth: np.ndarray) -> np.ndarray: + """Interpolate the horizon elevation at each solar azimuth. + + ``userhorizon`` follows the PVGIS convention: elevations in degrees at equally + spaced azimuths clockwise from north, the first entry being due north. The + profile wraps around, so the first entry is repeated at 360 degrees. + """ + horizon = np.asarray(userhorizon, dtype=float) + count = len(horizon) + azimuths = np.arange(count, dtype=float) * (360.0 / count) + return np.interp( + np.asarray(solar_azimuth, dtype=float) % 360.0, + np.append(azimuths, 360.0), + np.append(horizon, horizon[0]), + ) + + @staticmethod + def _tracked_orientation(plane: Any, solpos: pd.DataFrame) -> tuple[Any, Any]: + """Return (surface_tilt, surface_azimuth) honouring the plane's tracking type.""" + tilt = float(plane.surface_tilt if plane.surface_tilt is not None else 30.0) + azimuth = float(plane.surface_azimuth if plane.surface_azimuth is not None else 180.0) + tracking = plane.trackingtype + + if tracking in (None, 0): + return tilt, azimuth + + apparent_zenith = solpos["apparent_zenith"] + solar_azimuth = solpos["azimuth"] + + if tracking == 2: + # Two-axis: the plane always faces the sun. Below the horizon the angles are + # meaningless, so park the plane flat and let the zero irradiance do the rest. + return apparent_zenith.clip(lower=0.0, upper=90.0), solar_azimuth + if tracking == 3: + # Vertical axis: fixed tilt, azimuth follows the sun. + return tilt, solar_azimuth + if tracking in (1, 4, 5): + # Horizontal N-S (1), horizontal E-W (4), inclined N-S (5). + axis_tilt = tilt if tracking == 5 else 0.0 + axis_azimuth = 90.0 if tracking == 4 else 0.0 + tracker = pvlib.tracking.singleaxis( + apparent_zenith=apparent_zenith, + solar_azimuth=solar_azimuth, + axis_tilt=axis_tilt, + axis_azimuth=axis_azimuth, + max_angle=90, + backtrack=False, + ) + return ( + tracker["surface_tilt"].fillna(axis_tilt), + tracker["surface_azimuth"].fillna(axis_azimuth), + ) + + logger.warning( + f"Unsupported trackingtype {tracking} for local pvforecast, treating plane as fixed." + ) + return tilt, azimuth + + # ------------------------------------------------------------------ pv model + + def _plane_power( + self, + plane: Any, + weather: pd.DataFrame, + solpos: pd.DataFrame, + dni_extra: pd.Series, + airmass: pd.Series, + ) -> tuple[pd.Series, pd.Series]: + """Run the pvlib chain for one plane, returning (dc_power_w, ac_power_w).""" + settings = self._settings + + peakpower_kw = plane.peakpower + if peakpower_kw is None: + logger.warning("Plane without peakpower skipped by local pvforecast.") + zero = pd.Series(0.0, index=weather.index) + return zero, zero + pdc0 = float(peakpower_kw) * 1000.0 + + ghi = weather["shortwave_radiation"] + dhi = weather["diffuse_radiation"] + dni = weather["direct_normal_irradiance"] + + # Horizon shading kills the beam component; what is left of the global is the + # diffuse. Do this before transposition so the sky model sees consistent inputs. + if plane.userhorizon: + horizon = self._horizon_elevation(plane.userhorizon, solpos["azimuth"].to_numpy()) + shaded = solpos["apparent_elevation"].to_numpy() < horizon + dni = dni.where(~shaded, 0.0) + ghi = ghi.where(~shaded, dhi) + + surface_tilt, surface_azimuth = self._tracked_orientation(plane, solpos) + albedo = plane.albedo if plane.albedo is not None else settings.albedo + + poa = pvlib.irradiance.get_total_irradiance( + surface_tilt=surface_tilt, + surface_azimuth=surface_azimuth, + solar_zenith=solpos["apparent_zenith"], + solar_azimuth=solpos["azimuth"], + dni=dni, + ghi=ghi, + dhi=dhi, + dni_extra=dni_extra, + airmass=airmass, + albedo=float(albedo), + model=settings.transposition_model, + ) + poa_global = poa["poa_global"].fillna(0.0).clip(lower=0.0) + poa_direct = poa["poa_direct"].fillna(0.0).clip(lower=0.0) + poa_diffuse = poa["poa_diffuse"].fillna(0.0).clip(lower=0.0) + + if settings.apply_iam: + aoi = pvlib.irradiance.aoi( + surface_tilt, surface_azimuth, solpos["apparent_zenith"], solpos["azimuth"] + ) + iam = pvlib.iam.ashrae(aoi).fillna(0.0) + effective_irradiance = poa_direct * iam + poa_diffuse + else: + effective_irradiance = poa_global + + mounting = plane.mountingplace or "free" + temperature_params = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS["sapm"][ + MOUNTING_TEMPERATURE_MODEL.get(mounting, "open_rack_glass_glass") + ] + temp_cell = pvlib.temperature.sapm_cell( + poa_global=poa_global, + temp_air=weather["temperature_2m"], + wind_speed=weather["wind_speed_10m"], + **temperature_params, + ) + + dc_power = pvlib.pvsystem.pvwatts_dc( + effective_irradiance=effective_irradiance, + temp_cell=temp_cell, + pdc0=pdc0, + gamma_pdc=settings.temperature_coefficient / 100.0, + ) + # `loss` is the PVGIS-style lump of soiling, mismatch, wiring and ageing. + loss = plane.loss if plane.loss is not None else 0.0 + dc_power = (dc_power * (1.0 - float(loss) / 100.0)).fillna(0.0).clip(lower=0.0) + + paco = plane.inverter_paco + eta = settings.inverter_efficiency + if paco is None: + # No inverter rating configured: apply efficiency but do not clip. + ac_power = dc_power * eta + else: + ac_power = pd.Series( + pvlib.inverter.pvwatts( + pdc=dc_power.to_numpy(), + pdc0=float(paco) / eta, + eta_inv_nom=eta, + ), + index=dc_power.index, + ) + ac_power = ac_power.fillna(0.0).clip(lower=0.0) + + return dc_power, ac_power + + def _forecast_frame(self, data: Any, calibrate: bool = True) -> pd.DataFrame: + """Run the full chain and return a frame with dc/ac power indexed as EOS records. + + Args: + data: The Open-Meteo response. + calibrate: Apply the measurement-fitted correction when it is enabled and + fittable. Pass False to obtain the raw model output, which is what the + fit itself compares against. + """ + settings = self._settings + weather = self._weather_frame(data) + if weather.empty: + return pd.DataFrame(columns=["dc_power", "ac_power"]) + + location = pvlib.location.Location( + latitude=float(self.config.general.latitude), + longitude=float(self.config.general.longitude), + tz="UTC", + altitude=data.get("elevation"), + ) + + # Open-Meteo stamps an interval mean with the interval end, so the sun position + # that produced it sits half an interval earlier. + interval = pd.Timedelta(minutes=settings.resolution_minutes) + solar_times = weather.index - interval / 2 + + solpos_solar = location.get_solarposition(solar_times) + solpos = solpos_solar.set_axis(weather.index) + dni_extra = pd.Series( + np.asarray(pvlib.irradiance.get_extra_radiation(solar_times), dtype=float), + index=weather.index, + ) + airmass = pd.Series( + location.get_airmass(solar_times, solar_position=solpos_solar)[ + "airmass_relative" + ].to_numpy(), + index=weather.index, + ) + + total_dc = pd.Series(0.0, index=weather.index) + total_ac = pd.Series(0.0, index=weather.index) + for plane in self.config.pvforecast.planes or []: + dc_power, ac_power = self._plane_power(plane, weather, solpos, dni_extra, airmass) + total_dc = total_dc.add(dc_power, fill_value=0.0) + total_ac = total_ac.add(ac_power, fill_value=0.0) + + frame = pd.DataFrame( + { + "dc_power": total_dc, + "ac_power": total_ac, + "solar_elevation": solpos["apparent_elevation"].to_numpy(), + "solar_azimuth": solpos["azimuth"].to_numpy(), + } + ) + # Relabel from Open-Meteo's interval-end stamps to the interval-start stamps + # that EOS records use. + if settings.shift_to_interval_start: + frame.index = frame.index - interval + + if calibrate: + # The fit needs the raw model to compare against, which is exactly `frame`. + calibration = self._fit_calibration(frame) + if calibration is not None: + _, factors = calibration + frame = self._apply_calibration(frame, factors, self._installed_ac_capacity_w()) + return frame + + # ------------------------------------------------------------ calibration + + def _installed_ac_capacity_w(self) -> float: + """Rough installed AC capacity, used only to threshold near-zero intervals.""" + total = 0.0 + for plane in self.config.pvforecast.planes or []: + if plane.inverter_paco is not None: + total += float(plane.inverter_paco) + elif plane.peakpower is not None: + total += float(plane.peakpower) * 1000.0 + return total + + def _fit_calibration(self, frame: pd.DataFrame) -> Optional[tuple[float, np.ndarray]]: + """Fit correction factors from measured PV production against the model. + + The comparison runs on past intervals, where the Open-Meteo rows are analysed + rather than forecast weather. That is deliberate: it isolates the error of the + *PV model* (wrong kWp, soiling, degradation, shading the horizon profile misses) + from the error of the *weather forecast*, and only the former is systematic + enough to correct. + + Returns: + (global_factor, per_azimuth_bin_factors) or None when there is not enough + data to fit anything. + """ + settings = self._settings + if not settings.calibration_enabled: + return None + if not self.config.measurement.pv_production_emr_keys: + logger.info( + "PVForecastAkkudoktorLocal calibration is enabled but " + "measurement.pv_production_emr_keys is not configured - skipping." + ) + return None + + measurement = self.measurement + if measurement.max_datetime is None or measurement.min_datetime is None: + logger.info("PVForecastAkkudoktorLocal calibration: no PV measurements yet - skipping.") + return None + + interval = to_duration("1 hour") + end = measurement.max_datetime.start_of("hour") + start = end.subtract(days=settings.calibration_days) + if compare_datetimes(start, measurement.min_datetime).lt: + start = measurement.min_datetime.start_of("hour").add(hours=1) + # The model side only exists for the weather window that was requested. + model_start = to_datetime(frame.index[0].to_pydatetime()) + if compare_datetimes(start, model_start).lt: + start = model_start.start_of("hour").add(hours=1) + if compare_datetimes(start, end).ge: + logger.info("PVForecastAkkudoktorLocal calibration: measurement window too short - skipping.") + return None + + measured_kwh = np.asarray( + measurement.pv_production_total_kwh( + start_datetime=start, end_datetime=end, interval=interval + ), + dtype=float, + ) + if measured_kwh.size == 0 or not np.isfinite(measured_kwh).any(): + logger.info("PVForecastAkkudoktorLocal calibration: no usable PV measurements - skipping.") + return None + + # Model side on the same hourly grid. `ac_power` is a mean power per interval, + # so the hourly mean in W is directly the hourly energy in Wh. + hourly = frame[["ac_power", "solar_azimuth"]].resample("1h").mean() + grid = pd.date_range( + start=pd.Timestamp(start.in_timezone("UTC").isoformat()), + periods=len(measured_kwh), + freq="1h", + ) + hourly = hourly.reindex(grid) + modelled_kwh = hourly["ac_power"].to_numpy(dtype=float) / 1000.0 + azimuth = hourly["solar_azimuth"].to_numpy(dtype=float) + + # Only fit where the model says something meaningful is being produced. Dawn and + # dusk intervals otherwise dominate the ratio with noise. + floor_kwh = max(0.02 * self._installed_ac_capacity_w() / 1000.0, 0.05) + usable = ( + np.isfinite(modelled_kwh) + & np.isfinite(measured_kwh) + & np.isfinite(azimuth) + & (modelled_kwh > floor_kwh) + & (measured_kwh >= 0.0) + ) + if usable.sum() < 12: + logger.info( + f"PVForecastAkkudoktorLocal calibration: only {int(usable.sum())} usable hours - skipping." + ) + return None + + modelled_kwh = modelled_kwh[usable] + measured_kwh = measured_kwh[usable] + azimuth = azimuth[usable] + + model_total = float(modelled_kwh.sum()) + if model_total <= 0.0: + return None + global_factor = float( + np.clip( + measured_kwh.sum() / model_total, + settings.calibration_min_factor, + settings.calibration_max_factor, + ) + ) + + bin_degrees = settings.calibration_azimuth_bin_degrees + if bin_degrees <= 0: + logger.info( + f"PVForecastAkkudoktorLocal calibration: global factor {global_factor:.3f} " + f"from {usable.sum()} hours." + ) + return global_factor, np.array([global_factor]) + + bin_count = max(1, int(round(360 / bin_degrees))) + bin_index = np.clip((azimuth % 360.0) / (360.0 / bin_count), 0, bin_count - 1).astype(int) + + # Weight each bin by its modelled energy and shrink toward the global factor, so + # a thinly sampled bin cannot swing the forecast on its own. + prior = settings.calibration_prior_kwh + factors = np.full(bin_count, global_factor, dtype=float) + for b in range(bin_count): + in_bin = bin_index == b + weight = float(modelled_kwh[in_bin].sum()) + if weight <= 0.0: + continue + raw = float(measured_kwh[in_bin].sum()) / weight + factors[b] = np.clip( + (weight * raw + prior * global_factor) / (weight + prior), + settings.calibration_min_factor, + settings.calibration_max_factor, + ) + + # Report how much of the bias the fit actually removes on its own training window. + before = float(np.abs(modelled_kwh - measured_kwh).mean()) + after = float(np.abs(modelled_kwh * factors[bin_index] - measured_kwh).mean()) + logger.info( + f"PVForecastAkkudoktorLocal calibration over {usable.sum()} h: global factor " + f"{global_factor:.3f}, {bin_count} azimuth bins, " + f"MAE {before:.3f} -> {after:.3f} kWh/h" + ) + return global_factor, factors + + @staticmethod + def _apply_calibration( + frame: pd.DataFrame, factors: np.ndarray, ac_cap_w: float + ) -> pd.DataFrame: + """Scale the modelled power by the per-azimuth factor of each interval.""" + bin_count = len(factors) + azimuth = frame["solar_azimuth"].to_numpy(dtype=float) + bin_index = np.clip( + np.nan_to_num(azimuth % 360.0) / (360.0 / bin_count), 0, bin_count - 1 + ).astype(int) + scale = factors[bin_index] + frame = frame.copy() + frame["dc_power"] = frame["dc_power"] * scale + frame["ac_power"] = frame["ac_power"] * scale + if ac_cap_w > 0.0: + # A factor above 1 must not push the plant past its inverters. + frame["ac_power"] = frame["ac_power"].clip(upper=ac_cap_w) + return frame + + # ------------------------------------------------------------------ update + + def _update_data(self, force_update: Optional[bool] = False) -> None: + """Compute the PV forecast and store it as PVForecastDataRecord entries.""" + if not self.enabled(): + logger.info("PVForecastAkkudoktorLocal is disabled, skipping update.") + return + + if not self.config.pvforecast.planes: + error_msg = "Requested PV forecast, but no planes configured." + logger.error(f"Configuration error: {error_msg}") + raise ValueError(error_msg) + + data = self._request_forecast(force_update=force_update) # type: ignore[call-arg] + frame = self._forecast_frame(data) + if frame.empty: + logger.warning("Open-Meteo returned no weather rows for local pvforecast.") + return + + for timestamp, row in frame.iterrows(): + self.update_value( + to_datetime(timestamp.to_pydatetime()), + { + "pvforecast_dc_power": round(float(row["dc_power"]), 1), + "pvforecast_ac_power": round(float(row["ac_power"]), 1), + }, + ) + + logger.debug( + f"Updated local pvforecast: {len(frame)} records at " + f"{self._settings.resolution_minutes} min over " + f"{len(self.config.pvforecast.planes)} plane(s)." + ) + self.update_datetime = to_datetime(in_timezone=self.config.general.timezone) + + +# Example usage +if __name__ == "__main__": + pv = PVForecastAkkudoktorLocal() + pv._update_data() diff --git a/src/akkudoktoreos/prediction/weatheropenmeteo.py b/src/akkudoktoreos/prediction/weatheropenmeteo.py index e6957da3..23948e43 100644 --- a/src/akkudoktoreos/prediction/weatheropenmeteo.py +++ b/src/akkudoktoreos/prediction/weatheropenmeteo.py @@ -38,9 +38,9 @@ WeatherDataOpenMeteoMapping: List[Tuple[str, Optional[str], Optional[Union[str, ("wind_direction_10m", "Wind Direction (°)", 1), ("wind_gusts_10m", "Wind Gust Speed (kmph)", 3.6), # m/s to km/h ("shortwave_radiation", "Global Horizontal Irradiance (W/m2)", 1), - ("direct_radiation", "Direct Normal Irradiance (W/m2)", 1), + ("direct_radiation", None, None), # beam on the horizontal plane, not DNI ("diffuse_radiation", "Diffuse Horizontal Irradiance (W/m2)", 1), - ("direct_normal_irradiance", None, None), + ("direct_normal_irradiance", "Direct Normal Irradiance (W/m2)", 1), ("global_tilted_irradiance", None, None), ("terrestrial_radiation", None, None), ("shortwave_radiation_instant", None, None), @@ -148,7 +148,7 @@ class WeatherOpenMeteo(WeatherProvider): "wind_direction_10m", "wind_gusts_10m", "shortwave_radiation", # GHI - "direct_radiation", # DNI + "direct_normal_irradiance", # DNI "diffuse_radiation", # DHI "dew_point_2m", "apparent_temperature", @@ -157,6 +157,9 @@ class WeatherOpenMeteo(WeatherProvider): "sunshine_duration", ], "timezone": self.config.general.timezone, + # The mapping table converts m/s -> km/h, so ask for m/s explicitly instead + # of Open-Meteo's km/h default. + "wind_speed_unit": "ms", } # Calculate the number of days between start and end diff --git a/tests/test_prediction.py b/tests/test_prediction.py index 06177917..9d57a571 100644 --- a/tests/test_prediction.py +++ b/tests/test_prediction.py @@ -27,6 +27,7 @@ from akkudoktoreos.prediction.prediction import ( from akkudoktoreos.prediction.pvforecastakkudoktor import PVForecastAkkudoktor from akkudoktoreos.prediction.pvforecastforecastsolar import PVForecastForecastSolar from akkudoktoreos.prediction.pvforecastimport import PVForecastImport +from akkudoktoreos.prediction.pvforecastakkudoktorlocal import PVForecastAkkudoktorLocal from akkudoktoreos.prediction.pvforecastpvnode import PVForecastPVNode from akkudoktoreos.prediction.pvforecastsolcast import PVForecastSolcast from akkudoktoreos.prediction.pvforecastvrm import PVForecastVrm @@ -68,6 +69,7 @@ def forecast_providers(): PVForecastForecastSolar(), PVForecastSolcast(), PVForecastImport(), + PVForecastAkkudoktorLocal(), WeatherBrightSky(), WeatherClearOutside(), WeatherOpenMeteo(), @@ -126,10 +128,11 @@ def test_provider_sequence(prediction): assert isinstance(prediction.providers[19], PVForecastForecastSolar) assert isinstance(prediction.providers[20], PVForecastSolcast) assert isinstance(prediction.providers[21], PVForecastImport) - assert isinstance(prediction.providers[22], WeatherBrightSky) - assert isinstance(prediction.providers[23], WeatherClearOutside) - assert isinstance(prediction.providers[24], WeatherOpenMeteo) - assert isinstance(prediction.providers[25], WeatherImport) + assert isinstance(prediction.providers[22], PVForecastAkkudoktorLocal) + assert isinstance(prediction.providers[23], WeatherBrightSky) + assert isinstance(prediction.providers[24], WeatherClearOutside) + assert isinstance(prediction.providers[25], WeatherOpenMeteo) + assert isinstance(prediction.providers[26], WeatherImport) def test_provider_by_id(prediction, forecast_providers): diff --git a/tests/test_pvforecastakkudoktorlocal.py b/tests/test_pvforecastakkudoktorlocal.py new file mode 100644 index 00000000..151c7c69 --- /dev/null +++ b/tests/test_pvforecastakkudoktorlocal.py @@ -0,0 +1,291 @@ +"""Tests for the native (pvlib) PV forecast provider.""" + +from unittest.mock import patch + +import numpy as np +import pandas as pd +import pendulum +import pvlib +import pytest + +from akkudoktoreos.core.coreabc import get_measurement +from akkudoktoreos.prediction.pvforecastakkudoktorlocal import ( + PVForecastAkkudoktorLocal, + PVForecastAkkudoktorLocalCommonSettings, +) + +LATITUDE = 52.52 +LONGITUDE = 13.405 + +# A window that starts well before `START` so the calibration fit has past data. +WINDOW_START = pendulum.datetime(2025, 6, 1, 0, 0, tz="UTC") +WINDOW_END = pendulum.datetime(2025, 6, 20, 0, 0, tz="UTC") +START = pendulum.datetime(2025, 6, 15, 0, 0, tz="UTC") + + +def synthetic_openmeteo(resolution_minutes: int = 15, models: list[str] | None = None) -> dict: + """Build an Open-Meteo-shaped response from a pvlib clear-sky series. + + Open-Meteo stamps an interval mean with the interval END, and that is what the + provider expects, so the values are generated at those stamps directly. + """ + freq = f"{resolution_minutes}min" + index = pd.date_range( + start=WINDOW_START.format("YYYY-MM-DD HH:mm"), + end=WINDOW_END.format("YYYY-MM-DD HH:mm"), + freq=freq, + tz="UTC", + ) + location = pvlib.location.Location(LATITUDE, LONGITUDE, tz="UTC", altitude=37.0) + clearsky = location.get_clearsky(index, model="ineichen") + + block = "minutely_15" if resolution_minutes == 15 else "hourly" + values = { + "shortwave_radiation": clearsky["ghi"].round(1).tolist(), + "diffuse_radiation": clearsky["dhi"].round(1).tolist(), + "direct_normal_irradiance": clearsky["dni"].round(1).tolist(), + "temperature_2m": [20.0] * len(index), + "relative_humidity_2m": [50.0] * len(index), + "wind_speed_10m": [2.0] * len(index), + } + + data: dict = {"elevation": 37.0} + payload = {"time": [t.strftime("%Y-%m-%dT%H:%M") for t in index]} + if models: + # Multi-model requests come back with one suffixed series per member. + for name in models: + for key, series in values.items(): + payload[f"{key}_{name}"] = series + else: + payload.update(values) + data[block] = payload + return data + + +@pytest.fixture +def pvforecast_instance(config_eos): + config_eos.merge_settings_from_dict( + { + "general": {"latitude": LATITUDE, "longitude": LONGITUDE, "timezone": "UTC"}, + "prediction": {"hours": 96, "historic_hours": 48}, + "pvforecast": { + "provider": "PVForecastAkkudoktorLocal", + "planes": [ + { + "surface_tilt": 30.0, + "surface_azimuth": 180.0, + "peakpower": 10.0, + "inverter_paco": 10000, + "loss": 14.0, + } + ], + "provider_settings": {"PVForecastAkkudoktorLocal": {"resolution_minutes": 15}}, + }, + } + ) + return PVForecastAkkudoktorLocal(config=config_eos.load, start_datetime=START) + + +def test_provider_id(pvforecast_instance): + assert PVForecastAkkudoktorLocal.provider_id() == "PVForecastAkkudoktorLocal" + assert pvforecast_instance.enabled() is True + + +@pytest.mark.parametrize("value", [0, 5, 30, 61]) +def test_resolution_must_be_15_or_60(value): + with pytest.raises(ValueError, match="resolution_minutes must be 15 or 60"): + PVForecastAkkudoktorLocalCommonSettings(resolution_minutes=value) + + +def test_invalid_transposition_model(): + with pytest.raises(ValueError, match="Invalid transposition_model"): + PVForecastAkkudoktorLocalCommonSettings(transposition_model="nonsense") + + +def test_forecast_frame_is_quarter_hourly_and_plausible(pvforecast_instance): + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo()) + + assert not frame.empty + deltas = frame.index.to_series().diff().dropna().unique() + assert list(deltas) == [pd.Timedelta(minutes=15)] + + # A 10 kWp south-facing roof under clear June skies: below the inverter cap, + # but a substantial fraction of it. + peak = frame["ac_power"].max() + assert 5000.0 < peak <= 10000.0 + assert (frame["ac_power"] >= 0.0).all() + assert (frame["ac_power"] <= frame["dc_power"] + 1e-6).all() + + # Nights are dark. + midnight = frame.between_time("00:00", "01:00")["ac_power"] + assert midnight.max() == pytest.approx(0.0) + + +def test_records_are_shifted_to_interval_start(pvforecast_instance): + """Open-Meteo labels an interval by its end; EOS labels it by its start.""" + shifted = pvforecast_instance._forecast_frame(synthetic_openmeteo()) + + pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( + PVForecastAkkudoktorLocalCommonSettings(resolution_minutes=15, shift_to_interval_start=False) + ) + raw = pvforecast_instance._forecast_frame(synthetic_openmeteo()) + + assert raw.index[0] - shifted.index[0] == pd.Timedelta(minutes=15) + assert raw["ac_power"].to_numpy() == pytest.approx(shifted["ac_power"].to_numpy()) + + +def test_ensemble_members_are_averaged(pvforecast_instance): + """Several models in one request must be averaged, not dropped or duplicated.""" + single = pvforecast_instance._forecast_frame(synthetic_openmeteo()) + ensemble = pvforecast_instance._forecast_frame( + synthetic_openmeteo(models=["icon_seamless", "gfs_seamless", "ecmwf_ifs025"]) + ) + + # The synthetic members are identical, so the mean must reproduce the single run. + assert ensemble["ac_power"].to_numpy() == pytest.approx(single["ac_power"].to_numpy()) + + +def test_horizon_elevation_wraps_around_north(): + horizon = PVForecastAkkudoktorLocal._horizon_elevation( + [0.0, 10.0, 20.0, 30.0], np.array([0.0, 90.0, 180.0, 270.0, 359.999]) + ) + assert horizon[:4] == pytest.approx([0.0, 10.0, 20.0, 30.0]) + # Wrapping back to due north interpolates from 30 deg towards 0 deg. + assert horizon[4] == pytest.approx(0.0, abs=0.01) + + +def test_horizon_shading_reduces_yield(pvforecast_instance): + baseline = pvforecast_instance._forecast_frame(synthetic_openmeteo()) + + # A 40 deg wall all around blocks the beam for most of the day. + pvforecast_instance.config.pvforecast.planes[0].userhorizon = [40.0] * 12 + shaded = pvforecast_instance._forecast_frame(synthetic_openmeteo()) + + assert shaded["ac_power"].sum() < baseline["ac_power"].sum() * 0.9 + assert shaded["ac_power"].min() >= 0.0 + + # The low morning sun (below 40 deg elevation until ~07:00 UTC in June) is behind + # the wall, so its beam is gone entirely and only diffuse is left. + assert ( + shaded["ac_power"].between_time("05:00", "07:00").sum() + < baseline["ac_power"].between_time("05:00", "07:00").sum() * 0.5 + ) + + +def test_update_data_writes_records(pvforecast_instance): + with patch.object(PVForecastAkkudoktorLocal, "_request_forecast", return_value=synthetic_openmeteo()): + pvforecast_instance._update_data(force_update=True) + + assert len(pvforecast_instance.records) > 0 + record = pvforecast_instance.records[0] + assert record.pvforecast_ac_power is not None + assert record.pvforecast_dc_power is not None + + +def _feed_measurements( + instance: PVForecastAkkudoktorLocal, frame: pd.DataFrame, bias: float, key: str +) -> None: + """Write cumulative PV meter readings that are `bias` times the modelled power. + + Each test passes its own `key`. `Measurement` is database-backed, so clearing the + in-memory record list would not remove readings another test already stored. + """ + instance.config.measurement.pv_production_emr_keys = [key] + measurement = get_measurement() + + hourly = frame["ac_power"].resample("1h").mean() + hourly = hourly.loc[hourly.index < START] + + cumulative = 0.0 + for timestamp, power_w in hourly.items(): + measurement.update_value( + pendulum.instance(timestamp.to_pydatetime()), key, round(cumulative, 6) + ) + cumulative += float(power_w) * bias / 1000.0 + # Closing reading so the last interval has a difference to work with. + measurement.update_value( + pendulum.instance(hourly.index[-1].to_pydatetime()).add(hours=1), + key, + round(cumulative, 6), + ) + + +def test_calibration_is_off_by_default(pvforecast_instance): + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + assert pvforecast_instance._fit_calibration(frame) is None + + +def test_calibration_skips_without_measurement_keys(pvforecast_instance): + pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( + PVForecastAkkudoktorLocalCommonSettings(calibration_enabled=True) + ) + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + assert pvforecast_instance._fit_calibration(frame) is None + + +def test_calibration_recovers_a_systematic_bias(pvforecast_instance): + """A plant that consistently delivers 80% of the model must be corrected to 0.8.""" + pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( + PVForecastAkkudoktorLocalCommonSettings( + calibration_enabled=True, + calibration_days=14, + calibration_azimuth_bin_degrees=0, + ) + ) + + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + _feed_measurements(pvforecast_instance, frame, bias=0.8, key="pv_bias_emr") + + calibration = pvforecast_instance._fit_calibration(frame) + assert calibration is not None + global_factor, factors = calibration + assert global_factor == pytest.approx(0.8, abs=0.03) + assert factors == pytest.approx([global_factor]) + + corrected = pvforecast_instance._apply_calibration(frame, factors, 10000.0) + assert corrected["ac_power"].sum() == pytest.approx(frame["ac_power"].sum() * global_factor) + + +def test_calibration_factor_is_clamped(pvforecast_instance): + """A wildly wrong meter must not be allowed to swing the forecast.""" + pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( + PVForecastAkkudoktorLocalCommonSettings( + calibration_enabled=True, + calibration_days=14, + calibration_azimuth_bin_degrees=0, + calibration_min_factor=0.9, + calibration_max_factor=1.1, + ) + ) + + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + _feed_measurements(pvforecast_instance, frame, bias=0.2, key="pv_clamp_emr") + + calibration = pvforecast_instance._fit_calibration(frame) + assert calibration is not None + global_factor, _ = calibration + assert global_factor == pytest.approx(0.9) + + +def test_calibration_respects_the_inverter_cap(pvforecast_instance): + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + corrected = PVForecastAkkudoktorLocal._apply_calibration(frame, np.array([1.5]), 10000.0) + assert corrected["ac_power"].max() <= 10000.0 + 1e-6 + + +def test_forecast_frame_applies_the_calibration(pvforecast_instance): + """The correction must reach every caller of the chain, not just `_update_data`.""" + pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( + PVForecastAkkudoktorLocalCommonSettings( + calibration_enabled=True, + calibration_days=14, + calibration_azimuth_bin_degrees=0, + ) + ) + + raw = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + _feed_measurements(pvforecast_instance, raw, bias=0.8, key="pv_chain_emr") + + calibrated = pvforecast_instance._forecast_frame(synthetic_openmeteo()) + ratio = calibrated["ac_power"].sum() / raw["ac_power"].sum() + assert ratio == pytest.approx(0.8, abs=0.03) diff --git a/tests/test_weatheropenmeteo.py b/tests/test_weatheropenmeteo.py index 7e0136a4..6245754d 100644 --- a/tests/test_weatheropenmeteo.py +++ b/tests/test_weatheropenmeteo.py @@ -131,7 +131,7 @@ def test_request_forecast(mock_get, provider, sample_openmeteo_1_json): assert "time" in openmeteo_data["hourly"] assert "temperature_2m" in openmeteo_data["hourly"] assert "shortwave_radiation" in openmeteo_data["hourly"] # GHI - assert "direct_radiation" in openmeteo_data["hourly"] # DNI + assert "direct_normal_irradiance" in openmeteo_data["hourly"] # DNI assert "diffuse_radiation" in openmeteo_data["hourly"] # DHI @@ -161,7 +161,7 @@ def test_update_data(mock_get, provider, sample_openmeteo_1_json, cache_store): # Get the first record and check for irradiance values value_datetime = to_datetime("2026-03-04 09:00:00+01:00", in_timezone="Europe/Berlin") assert provider.key_to_value("weather_ghi", target_datetime=start_datetime) == 21.8 - assert provider.key_to_value("weather_dni", target_datetime=start_datetime) == 1.2 + assert provider.key_to_value("weather_dni", target_datetime=start_datetime) == 17.9 assert provider.key_to_value("weather_dhi", target_datetime=start_datetime) == 20.5 @@ -176,18 +176,22 @@ def test_openmeteo_radiation_mapping(provider): from akkudoktoreos.prediction.weatheropenmeteo import WeatherDataOpenMeteoMapping radiation_keys = [item[0] for item in WeatherDataOpenMeteoMapping - if item[0] in ['shortwave_radiation', 'direct_radiation', 'diffuse_radiation']] + if item[0] in ['shortwave_radiation', 'direct_normal_irradiance', + 'diffuse_radiation']] assert 'shortwave_radiation' in radiation_keys - assert 'direct_radiation' in radiation_keys + assert 'direct_normal_irradiance' in radiation_keys assert 'diffuse_radiation' in radiation_keys - # Verify they map to correct descriptions + # Verify they map to correct descriptions. Open-Meteo's `direct_radiation` is beam + # irradiance on the HORIZONTAL plane, so it must not be mapped to DNI. for key, desc, _ in WeatherDataOpenMeteoMapping: if key == 'shortwave_radiation': assert desc == "Global Horizontal Irradiance (W/m2)" - elif key == 'direct_radiation': + elif key == 'direct_normal_irradiance': assert desc == "Direct Normal Irradiance (W/m2)" + elif key == 'direct_radiation': + assert desc is None elif key == 'diffuse_radiation': assert desc == "Diffuse Horizontal Irradiance (W/m2)" diff --git a/tests/testdata/weatherforecast_openmeteo_1.json b/tests/testdata/weatherforecast_openmeteo_1.json index 648fa307..0fb57486 100644 --- a/tests/testdata/weatherforecast_openmeteo_1.json +++ b/tests/testdata/weatherforecast_openmeteo_1.json @@ -8,7 +8,7 @@ "elevation": 291.0, "hourly_units": { "time": "iso8601", - "temperature_2m": "\u00b0C", + "temperature_2m": "°C", "relative_humidity_2m": "%", "precipitation": "mm", "rain": "mm", @@ -19,16 +19,17 @@ "pressure_msl": "hPa", "surface_pressure": "hPa", "wind_speed_10m": "km/h", - "wind_direction_10m": "\u00b0", + "wind_direction_10m": "°", "wind_gusts_10m": "km/h", - "shortwave_radiation": "W/m\u00b2", - "direct_radiation": "W/m\u00b2", - "diffuse_radiation": "W/m\u00b2", - "dew_point_2m": "\u00b0C", - "apparent_temperature": "\u00b0C", + "shortwave_radiation": "W/m²", + "direct_radiation": "W/m²", + "diffuse_radiation": "W/m²", + "dew_point_2m": "°C", + "apparent_temperature": "°C", "precipitation_probability": "%", "visibility": "m", - "sunshine_duration": "s" + "sunshine_duration": "s", + "direct_normal_irradiance": "W/m²" }, "hourly": { "time": [ @@ -1658,6 +1659,80 @@ 0.0, 0.0, 0.0 + ], + "direct_normal_irradiance": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 17.9, + 62.9, + 174.7, + 669.0, + 735.6, + 765.5, + 576.9, + 433.0, + 631.8, + 465.1, + 222.8, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 11.0, + 44.4, + 76.8, + 394.1, + 736.5, + 774.1, + 772.9, + 738.2, + 656.6, + 502.2, + 225.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 6.4, + 54.1, + 205.1, + 490.6, + 735.7, + 767.7, + 771.1, + 731.0, + 648.3, + 498.4, + 224.8, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 ] } -} \ No newline at end of file +} diff --git a/tests/testdata/weatherforecast_openmeteo_2.json b/tests/testdata/weatherforecast_openmeteo_2.json index 0b2b7ee4..68adf6c3 100644 --- a/tests/testdata/weatherforecast_openmeteo_2.json +++ b/tests/testdata/weatherforecast_openmeteo_2.json @@ -231,7 +231,7 @@ "weather_pressure": 10.26, "weather_ozone": null, "weather_ghi": 0.0, - "weather_dni": 0.0, + "weather_dni": 17.9, "weather_dhi": 0.0 }, { @@ -257,7 +257,7 @@ "weather_pressure": 10.265, "weather_ozone": null, "weather_ghi": 21.8, - "weather_dni": 1.2, + "weather_dni": 62.9, "weather_dhi": 20.5 }, { @@ -283,7 +283,7 @@ "weather_pressure": 10.269000000000002, "weather_ozone": null, "weather_ghi": 101.2, - "weather_dni": 13.8, + "weather_dni": 174.7, "weather_dhi": 87.5 }, { @@ -309,7 +309,7 @@ "weather_pressure": 10.262, "weather_ozone": null, "weather_ghi": 207.0, - "weather_dni": 61.5, + "weather_dni": 669.0, "weather_dhi": 145.5 }, { @@ -335,7 +335,7 @@ "weather_pressure": 10.257000000000001, "weather_ozone": null, "weather_ghi": 407.5, - "weather_dni": 304.2, + "weather_dni": 735.6, "weather_dhi": 103.2 }, { @@ -361,7 +361,7 @@ "weather_pressure": 10.252, "weather_ozone": null, "weather_ghi": 487.2, - "weather_dni": 382.5, + "weather_dni": 765.5, "weather_dhi": 104.8 }, { @@ -387,7 +387,7 @@ "weather_pressure": 10.247, "weather_ozone": null, "weather_ghi": 519.5, - "weather_dni": 416.0, + "weather_dni": 576.9, "weather_dhi": 103.5 }, { @@ -413,7 +413,7 @@ "weather_pressure": 10.235, "weather_ozone": null, "weather_ghi": 450.8, - "weather_dni": 302.0, + "weather_dni": 433.0, "weather_dhi": 148.8 }, { @@ -439,7 +439,7 @@ "weather_pressure": 10.23, "weather_ozone": null, "weather_ghi": 367.2, - "weather_dni": 199.8, + "weather_dni": 631.8, "weather_dhi": 167.5 }, { @@ -465,7 +465,7 @@ "weather_pressure": 10.228, "weather_ozone": null, "weather_ghi": 315.0, - "weather_dni": 228.5, + "weather_dni": 465.1, "weather_dhi": 86.5 }, { @@ -491,7 +491,7 @@ "weather_pressure": 10.227, "weather_ozone": null, "weather_ghi": 174.2, - "weather_dni": 107.5, + "weather_dni": 222.8, "weather_dhi": 66.8 }, { @@ -517,7 +517,7 @@ "weather_pressure": 10.23, "weather_ozone": null, "weather_ghi": 42.5, - "weather_dni": 17.8, + "weather_dni": 0.0, "weather_dhi": 24.8 }, { @@ -855,7 +855,7 @@ "weather_pressure": 10.278, "weather_ozone": null, "weather_ghi": 0.0, - "weather_dni": 0.0, + "weather_dni": 11.0, "weather_dhi": 0.0 }, { @@ -881,7 +881,7 @@ "weather_pressure": 10.287, "weather_ozone": null, "weather_ghi": 23.5, - "weather_dni": 0.8, + "weather_dni": 44.4, "weather_dhi": 22.8 }, { @@ -907,7 +907,7 @@ "weather_pressure": 10.288, "weather_ozone": null, "weather_ghi": 90.0, - "weather_dni": 10.0, + "weather_dni": 76.8, "weather_dhi": 80.0 }, { @@ -933,7 +933,7 @@ "weather_pressure": 10.286, "weather_ozone": null, "weather_ghi": 177.0, - "weather_dni": 27.5, + "weather_dni": 394.1, "weather_dhi": 149.5 }, { @@ -959,7 +959,7 @@ "weather_pressure": 10.279000000000002, "weather_ozone": null, "weather_ghi": 352.0, - "weather_dni": 181.5, + "weather_dni": 736.5, "weather_dhi": 170.5 }, { @@ -985,7 +985,7 @@ "weather_pressure": 10.272, "weather_ozone": null, "weather_ghi": 496.0, - "weather_dni": 387.2, + "weather_dni": 774.1, "weather_dhi": 108.8 }, { @@ -1011,7 +1011,7 @@ "weather_pressure": 10.267999999999999, "weather_ozone": null, "weather_ghi": 530.5, - "weather_dni": 425.0, + "weather_dni": 772.9, "weather_dhi": 105.5 }, { @@ -1037,7 +1037,7 @@ "weather_pressure": 10.265, "weather_ozone": null, "weather_ghi": 509.0, - "weather_dni": 408.8, + "weather_dni": 738.2, "weather_dhi": 100.2 }, { @@ -1063,7 +1063,7 @@ "weather_pressure": 10.261, "weather_ozone": null, "weather_ghi": 437.2, - "weather_dni": 344.5, + "weather_dni": 656.6, "weather_dhi": 92.8 }, { @@ -1089,7 +1089,7 @@ "weather_pressure": 10.257000000000001, "weather_ozone": null, "weather_ghi": 321.2, - "weather_dni": 240.8, + "weather_dni": 502.2, "weather_dhi": 80.5 }, { @@ -1115,7 +1115,7 @@ "weather_pressure": 10.257000000000001, "weather_ozone": null, "weather_ghi": 179.0, - "weather_dni": 118.5, + "weather_dni": 225.0, "weather_dhi": 60.5 }, { @@ -1141,7 +1141,7 @@ "weather_pressure": 10.26, "weather_ozone": null, "weather_ghi": 44.2, - "weather_dni": 19.0, + "weather_dni": 0.0, "weather_dhi": 25.2 }, { @@ -1479,7 +1479,7 @@ "weather_pressure": 10.290000000000001, "weather_ozone": null, "weather_ghi": 0.0, - "weather_dni": 0.0, + "weather_dni": 6.4, "weather_dhi": 0.0 }, { @@ -1505,7 +1505,7 @@ "weather_pressure": 10.293, "weather_ozone": null, "weather_ghi": 24.2, - "weather_dni": 0.5, + "weather_dni": 54.1, "weather_dhi": 23.8 }, { @@ -1531,7 +1531,7 @@ "weather_pressure": 10.295, "weather_ozone": null, "weather_ghi": 98.0, - "weather_dni": 12.5, + "weather_dni": 205.1, "weather_dhi": 85.5 }, { @@ -1557,7 +1557,7 @@ "weather_pressure": 10.294, "weather_ozone": null, "weather_ghi": 218.8, - "weather_dni": 74.6, + "weather_dni": 490.6, "weather_dhi": 144.2 }, { @@ -1583,7 +1583,7 @@ "weather_pressure": 10.288, "weather_ozone": null, "weather_ghi": 380.7, - "weather_dni": 228.8, + "weather_dni": 735.7, "weather_dhi": 151.8 }, { @@ -1609,7 +1609,7 @@ "weather_pressure": 10.278, "weather_ozone": null, "weather_ghi": 496.1, - "weather_dni": 391.0, + "weather_dni": 767.7, "weather_dhi": 105.1 }, { @@ -1635,7 +1635,7 @@ "weather_pressure": 10.267999999999999, "weather_ozone": null, "weather_ghi": 528.0, - "weather_dni": 425.8, + "weather_dni": 771.1, "weather_dhi": 102.2 }, { @@ -1661,7 +1661,7 @@ "weather_pressure": 10.263, "weather_ozone": null, "weather_ghi": 508.0, - "weather_dni": 412.0, + "weather_dni": 731.0, "weather_dhi": 96.0 }, { @@ -1687,7 +1687,7 @@ "weather_pressure": 10.261, "weather_ozone": null, "weather_ghi": 435.0, - "weather_dni": 345.0, + "weather_dni": 648.3, "weather_dhi": 90.0 }, { @@ -1713,7 +1713,7 @@ "weather_pressure": 10.258, "weather_ozone": null, "weather_ghi": 320.0, - "weather_dni": 241.0, + "weather_dni": 498.4, "weather_dhi": 79.0 }, { @@ -1739,7 +1739,7 @@ "weather_pressure": 10.255, "weather_ozone": null, "weather_ghi": 180.0, - "weather_dni": 120.0, + "weather_dni": 224.8, "weather_dhi": 60.0 }, { @@ -1765,7 +1765,7 @@ "weather_pressure": 10.255, "weather_ozone": null, "weather_ghi": 45.0, - "weather_dni": 20.0, + "weather_dni": 0.0, "weather_dhi": 25.0 }, { @@ -1956,4 +1956,4 @@ "keep_datetime": "2026-02-28 00:00:00+01:00", "total_hours": 48, "keep_hours": 48 -} \ No newline at end of file +}