diff --git a/CHANGELOG.md b/CHANGELOG.md index 6881fe62..e12c3fd9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -118,6 +118,15 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). - 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. +- The local PV provider's calibration now excludes probable outage and curtailment days instead of + learning them as permanent model losses: `calibration_outage_filter_enabled` (default on), + `calibration_outage_threshold`, `calibration_reference_days` and `calibration_min_healthy_days` + estimate the healthy plant ratio and fall back to the most recent healthy days. Calibration also + uses native 15-minute meter readings when every configured PV meter supplies them, interpolates + azimuth factors smoothly between bin centres instead of stepping, and normalizes the fitted + shape per forecast day so it redistributes energy without changing that day's kWh correction. + The default `calibration_azimuth_bin_degrees` moves from 15 to 45, which is what a typical + calibration window actually supports. - Separate the control horizon from the battery lookahead. `optimization.horizon_hours` remains the only span that receives control commands; the new `optimization.tail_horizon_hours` (default 48 h) is a forecast lookahead that never produces a command. In `AUTO` terminal-value @@ -219,6 +228,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). provider could 404 the whole load prediction. It now defaults to a cache-aware update and accepts an optional `force_update` flag in the request body for callers that still want to force. +- The local PV provider derived its calibration window from the measurement store as a whole + instead of from the configured PV production meters. A load meter reaching further than the PV + meter placed the window where no PV reading exists, so calibration silently fell back to hourly + fitting or skipped itself entirely. The window now follows the PV meters. +- `Measurement.load()` silently discarded every stored record. It validated the file into a + temporary `Measurement`, but `Measurement` is a singleton, so the "temporary" instance was the + already initialized one and the parsed records were dropped. The records are now validated + individually and inserted directly. - A rejected configuration update no longer damages the running configuration. `merge_settings_from_dict` validated the merged candidate only while reinitializing the singleton, so an invalid update could leave EOS half-updated. The candidate is validated first. diff --git a/docs/_generated/configpvforecast.md b/docs/_generated/configpvforecast.md index f6d1811c..401dadcb 100644 --- a/docs/_generated/configpvforecast.md +++ b/docs/_generated/configpvforecast.md @@ -206,12 +206,16 @@ | ---- | ---- | --------- | ------- | ----------- | | 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_azimuth_bin_degrees | `int` | `rw` | `45` | 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_min_healthy_days | `int` | `rw` | `3` | Minimum number of healthy days used for a fit. Older healthy days from the reference window are added when the recent window contains fewer. | +| calibration_outage_filter_enabled | `bool` | `rw` | `True` | Exclude days whose measured production is far below the recent healthy plant level. This prevents inverter, battery and curtailment events from being learned as permanent PV model losses. | +| calibration_outage_threshold | `float` | `rw` | `0.55` | A day is treated as unavailable when its measured/modelled energy ratio is below this fraction of the robust healthy reference ratio. | | 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. | +| calibration_reference_days | `int` | `rw` | `30` | Lookback used to distinguish healthy production from outages or curtailment. If the calibration window contains too few healthy days, the most recent healthy days from this reference window are used. | | 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`. | @@ -247,7 +251,11 @@ "shift_to_interval_start": true, "calibration_enabled": true, "calibration_days": 30, - "calibration_azimuth_bin_degrees": 15, + "calibration_reference_days": 30, + "calibration_outage_filter_enabled": true, + "calibration_outage_threshold": 0.55, + "calibration_min_healthy_days": 3, + "calibration_azimuth_bin_degrees": 45, "calibration_prior_kwh": 5.0, "calibration_min_factor": 0.5, "calibration_max_factor": 1.5 @@ -432,9 +440,9 @@ | Name | Type | Read-Only | Default | Description | | ---- | ---- | --------- | ------- | ----------- | +| PVForecastAkkudoktorLocal | `Optional[akkudoktoreos.prediction.pvforecastakkudoktorlocal.PVForecastAkkudoktorLocalCommonSettings]` | `rw` | `None` | PVForecastAkkudoktorLocal settings | | 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 | diff --git a/docs/akkudoktoreos/prediction.md b/docs/akkudoktoreos/prediction.md index 7c0ef7e3..2ad48225 100644 --- a/docs/akkudoktoreos/prediction.md +++ b/docs/akkudoktoreos/prediction.md @@ -822,6 +822,25 @@ factor is clamped to `[calibration_min_factor, calibration_max_factor]` so a bro either. The fitted factors and the resulting change in mean absolute error are logged at INFO level on every update. +By default, calibration also rejects probable outage or curtailment days. It estimates the +healthy plant ratio from `calibration_reference_days`, excludes days below +`calibration_outage_threshold` of that reference, and falls back to the most recent +`calibration_min_healthy_days` when the normal calibration window contains an outage. This keeps +a battery or inverter failure that limits PV to local demand from becoming a permanent forecast +loss. Set `calibration_outage_filter_enabled` to false only when measured curtailed production, +rather than available PV potential, is the intended prediction target. + +The measurement cadence controls the detail that can be learned. Hourly cumulative meter +readings calibrate hourly energy while the native Open-Meteo/pvlib chain continues to supply the +15-minute shape. If every configured PV meter supplies genuine 15-minute readings, calibration +automatically uses those native slots as well. It never interpolates hourly counters into an +invented quarter-hour profile. Azimuth factors are interpolated smoothly between bin centres so +they do not introduce steps into the EMS input curve. The shape fit uses all healthy days in the +reference window, while the global factor still follows the shorter recent window. Finally, the +shape is normalized per forecast day: it redistributes the calibrated energy across the day's +15-minute slots without changing that day's global kWh correction (unless the physical inverter +limit clips a peak). + Calibration requires `measurement.pv_production_emr_keys` to be configured and fed with cumulative PV production meter readings in kWh: diff --git a/scripts/pvforecast_backtest.py b/scripts/pvforecast_backtest.py index 02f06aa0..4a0cd363 100644 --- a/scripts/pvforecast_backtest.py +++ b/scripts/pvforecast_backtest.py @@ -46,8 +46,11 @@ 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}), + ( + "best_match", + {"weather_models": ["best_match"], "calibration_enabled": False}, + ), + ("ensemble", {"weather_models": ENSEMBLE, "calibration_enabled": False}), ("ensemble + calibration", {"weather_models": ENSEMBLE, "calibration_enabled": True}), ( "ensemble + calibration (global only)", @@ -57,8 +60,18 @@ VARIANTS: list[tuple[str, dict[str, Any]]] = [ "calibration_azimuth_bin_degrees": 0, }, ), - ("ensemble, isotropic sky", {"weather_models": ENSEMBLE, "transposition_model": "isotropic"}), - ("ensemble, no IAM", {"weather_models": ENSEMBLE, "apply_iam": False}), + ( + "ensemble, isotropic sky", + { + "weather_models": ENSEMBLE, + "transposition_model": "isotropic", + "calibration_enabled": False, + }, + ), + ( + "ensemble, no IAM", + {"weather_models": ENSEMBLE, "apply_iam": False, "calibration_enabled": False}, + ), ] @@ -94,6 +107,8 @@ def main(days: int, tilt: Optional[float], azimuth: Optional[float]) -> int: singletons_init() config = get_config() measurement = get_measurement() + if measurement.max_datetime is None: + measurement.load() if not config.measurement.pv_production_emr_keys: print( @@ -175,6 +190,11 @@ def main(days: int, tilt: Optional[float], azimuth: Optional[float]) -> int: "on, so their advantage here is optimistic. Re-run with a longer --days to see\n" "how much of it survives." ) + print( + "Outage or curtailment periods are excluded from calibration but remain in these\n" + "scores, because the script cannot prove the plant's availability without an\n" + "explicit availability measurement." + ) best_label, best, hours = min(rows, key=lambda row: row[1]["mae"]) print( diff --git a/src/akkudoktoreos/measurement/measurement.py b/src/akkudoktoreos/measurement/measurement.py index 3c6d1328..cab4cd14 100644 --- a/src/akkudoktoreos/measurement/measurement.py +++ b/src/akkudoktoreos/measurement/measurement.py @@ -6,6 +6,7 @@ data records for measurements. The measurements can be added programmatically or imported from a file or JSON string. """ +import json from pathlib import Path from typing import Any, Optional @@ -359,14 +360,12 @@ class Measurement(SingletonMixin, DataImportMixin, DataSequence): if not measurement_file_path.exists(): return False try: - # Validate into a temporary instance - loaded = self.__class__.model_validate_json( - measurement_file_path.read_text(encoding="utf-8") - ) - - # Explicitly add data records to the existing singleton - for record in loaded.records: - self.insert_by_datetime(record) + # Do not validate the complete Measurement model here. Measurement is a + # singleton, so constructing a temporary instance returns the already + # initialized singleton and silently discards the serialized records. + payload = json.loads(measurement_file_path.read_text(encoding="utf-8")) + for record_data in payload.get("records", []): + self.insert_by_datetime(MeasurementDataRecord.model_validate(record_data)) except Exception as e: logger.exception("Cannot load measurements") return True diff --git a/src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py b/src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py index d31b2e55..61f3da12 100644 --- a/src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py +++ b/src/akkudoktoreos/prediction/pvforecastakkudoktorlocal.py @@ -203,8 +203,56 @@ class PVForecastAkkudoktorLocalCommonSettings(SettingsBaseModel): "examples": [30, 14], }, ) + calibration_reference_days: int = Field( + default=30, + ge=3, + le=MAX_PAST_DAYS, + json_schema_extra={ + "description": ( + "Lookback used to distinguish healthy production from outages or " + "curtailment. If the calibration window contains too few healthy days, " + "the most recent healthy days from this reference window are used." + ), + "examples": [30, 14], + }, + ) + calibration_outage_filter_enabled: bool = Field( + default=True, + json_schema_extra={ + "description": ( + "Exclude days whose measured production is far below the recent healthy " + "plant level. This prevents inverter, battery and curtailment events from " + "being learned as permanent PV model losses." + ), + "examples": [True], + }, + ) + calibration_outage_threshold: float = Field( + default=0.55, + gt=0.0, + lt=1.0, + json_schema_extra={ + "description": ( + "A day is treated as unavailable when its measured/modelled energy ratio " + "is below this fraction of the robust healthy reference ratio." + ), + "examples": [0.55, 0.5], + }, + ) + calibration_min_healthy_days: int = Field( + default=3, + ge=1, + le=31, + json_schema_extra={ + "description": ( + "Minimum number of healthy days used for a fit. Older healthy days from " + "the reference window are added when the recent window contains fewer." + ), + "examples": [3], + }, + ) calibration_azimuth_bin_degrees: int = Field( - default=15, + default=45, ge=0, le=180, json_schema_extra={ @@ -212,7 +260,7 @@ class PVForecastAkkudoktorLocalCommonSettings(SettingsBaseModel): "Width of the solar-azimuth bins for the correction. 0 fits a single " "global factor only." ), - "examples": [15, 30, 0], + "examples": [45, 30, 15, 0], }, ) calibration_prior_kwh: float = Field( @@ -295,7 +343,12 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): 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) + calibration_lookback = settings.calibration_days + if settings.calibration_outage_filter_enabled: + calibration_lookback = max( + calibration_lookback, settings.calibration_reference_days + ) + past_days = max(past_days, calibration_lookback) return forecast_days, min(MAX_PAST_DAYS, past_days) @@ -305,7 +358,9 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): 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") + raise ValueError( + "PVForecastAkkudoktorLocal needs general.latitude and general.longitude" + ) settings = self._settings block = "minutely_15" if settings.resolution_minutes == 15 else "hourly" @@ -602,8 +657,14 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): # 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()) + global_factor, factors = calibration + frame = self._apply_calibration( + frame, + factors, + self._installed_ac_capacity_w(), + global_factor=global_factor, + timezone=self.config.general.timezone, + ) return frame # ------------------------------------------------------------ calibration @@ -618,6 +679,127 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): total += float(plane.peakpower) * 1000.0 return total + def _pv_measurement_window(self) -> Optional[tuple[Any, Any]]: + """First and last timestamp that actually carries PV production readings. + + The measurement store holds every meter, not just the PV ones. Load meters + routinely reach further than the PV meter in both directions, so deriving + the calibration window from the store as a whole would place it where no + PV reading exists - which silently skips calibration or downgrades it to + hourly fitting. + """ + measurement = self.measurement + if measurement.min_datetime is None or measurement.max_datetime is None: + return None + earliest: Any = None + latest: Any = None + for key in self.config.measurement.pv_production_emr_keys or []: + dates, _ = measurement.key_to_lists( + key=key, + start_datetime=measurement.min_datetime, + end_datetime=measurement.max_datetime.add(minutes=1), + ) + if not dates: + continue + if earliest is None or compare_datetimes(dates[0], earliest).lt: + earliest = dates[0] + if latest is None or compare_datetimes(dates[-1], latest).gt: + latest = dates[-1] + if earliest is None or latest is None: + return None + return earliest, latest + + def _calibration_interval_minutes(self, start: Any, end: Any) -> int: + """Use native forecast slots only when every PV meter resolves them. + + Interpolating an hourly cumulative meter onto quarter hours would create a + perfectly flat, but invented, intrahour profile. Fall back to hourly fitting + until all configured production meters actually provide native slot readings. + """ + native_minutes = self._settings.resolution_minutes + meter_resolutions: list[float] = [] + for key in self.config.measurement.pv_production_emr_keys or []: + dates, _ = self.measurement.key_to_lists( + key=key, start_datetime=start, end_datetime=end + ) + if len(dates) < 3: + return 60 + deltas = np.asarray( + [ + (dates[index] - dates[index - 1]).total_seconds() / 60.0 + for index in range(1, len(dates)) + if dates[index] > dates[index - 1] + ], + dtype=float, + ) + if deltas.size == 0: + return 60 + meter_resolutions.append(float(np.median(deltas))) + + if not meter_resolutions: + return 60 + if max(meter_resolutions) <= native_minutes * 1.5: + return native_minutes + return 60 + + def _fit_azimuth_factors( + self, + modelled_kwh: np.ndarray, + measured_kwh: np.ndarray, + azimuth: np.ndarray, + global_factor: float, + ) -> np.ndarray: + """Fit an energy-weighted intraday shape while preserving global energy.""" + settings = self._settings + bin_degrees = settings.calibration_azimuth_bin_degrees + if bin_degrees <= 0: + return 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) + + # Fit the shape as a residual around the independently determined global + # energy factor. Day-level availability filtering has already removed outages; + # energy sums now give productive intervals the influence relevant to the EMS. + # The prior shrinks sparse bins back toward a neutral relative factor of one. + relative_shape = np.ones(bin_count, dtype=float) + prior = settings.calibration_prior_kwh + for bin_number in range(bin_count): + in_bin = bin_index == bin_number + weight = float(modelled_kwh[in_bin].sum()) + if weight <= 0.0: + continue + measured_sum = float(measured_kwh[in_bin].sum()) + raw_shape = measured_sum / (weight * global_factor) + relative_shape[bin_number] = (weight * raw_shape + prior) / (weight + prior) + + factors = np.clip( + global_factor * relative_shape, + settings.calibration_min_factor, + settings.calibration_max_factor, + ) + + # Smooth interpolation changes the exact weighted mean of the bin-centre + # values. Renormalize after interpolation so shape correction cannot silently + # change the global kWh calibration. Re-clipping is iterated to respect bounds. + target_energy = global_factor * float(modelled_kwh.sum()) + for _ in range(8): + scale = self._interpolate_azimuth_factors(azimuth, factors) + corrected_energy = float(np.dot(modelled_kwh, scale)) + if corrected_energy <= 0.0: + break + correction = target_energy / corrected_energy + updated = np.clip( + factors * correction, + settings.calibration_min_factor, + settings.calibration_max_factor, + ) + if np.allclose(updated, factors, rtol=1e-6, atol=1e-8): + factors = updated + break + factors = updated + return factors + def _fit_calibration(self, frame: pd.DataFrame) -> Optional[tuple[float, np.ndarray]]: """Fit correction factors from measured PV production against the model. @@ -642,21 +824,35 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): return None measurement = self.measurement - if measurement.max_datetime is None or measurement.min_datetime is None: + pv_window = self._pv_measurement_window() + if pv_window is None: logger.info("PVForecastAkkudoktorLocal calibration: no PV measurements yet - skipping.") return None + pv_min_datetime, reference_end = pv_window - 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) + lookback_days = settings.calibration_days + if settings.calibration_outage_filter_enabled: + lookback_days = max(lookback_days, settings.calibration_reference_days) + reference_start = reference_end.subtract(days=lookback_days) + interval_minutes = self._calibration_interval_minutes(reference_start, reference_end) + interval_hours = interval_minutes / 60.0 + interval = to_duration(f"{interval_minutes} minutes") + end = reference_end.start_of("hour") + if interval_minutes < 60: + end = reference_end.start_of("minute").set( + minute=(reference_end.minute // interval_minutes) * interval_minutes + ) + start = end.subtract(days=lookback_days) + if compare_datetimes(start, pv_min_datetime).lt: + start = pv_min_datetime.start_of("minute").add(minutes=interval_minutes) # 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) + start = model_start.start_of("minute").add(minutes=interval_minutes) if compare_datetimes(start, end).ge: - logger.info("PVForecastAkkudoktorLocal calibration: measurement window too short - skipping.") + logger.info( + "PVForecastAkkudoktorLocal calibration: measurement window too short - skipping." + ) return None measured_kwh = np.asarray( @@ -666,24 +862,32 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): dtype=float, ) if measured_kwh.size == 0 or not np.isfinite(measured_kwh).any(): - logger.info("PVForecastAkkudoktorLocal calibration: no usable PV measurements - skipping.") + 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() + # Model side on the same grid as the real meter. Mean power is converted to + # interval energy below; hourly meters stay hourly and native 15-minute meters + # retain the shape that matters to the EMS. + samples_frame = ( + frame[["ac_power", "solar_azimuth"]].resample(f"{interval_minutes}min").mean() + ) grid = pd.date_range( start=pd.Timestamp(start.in_timezone("UTC").isoformat()), periods=len(measured_kwh), - freq="1h", + freq=f"{interval_minutes}min", ) - hourly = hourly.reindex(grid) - modelled_kwh = hourly["ac_power"].to_numpy(dtype=float) / 1000.0 - azimuth = hourly["solar_azimuth"].to_numpy(dtype=float) + samples_frame = samples_frame.reindex(grid) + modelled_kwh = samples_frame["ac_power"].to_numpy(dtype=float) / 1000.0 * interval_hours + azimuth = samples_frame["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) + floor_kwh = max( + 0.02 * self._installed_ac_capacity_w() / 1000.0 * interval_hours, + 0.05 * interval_hours, + ) usable = ( np.isfinite(modelled_kwh) & np.isfinite(measured_kwh) @@ -691,12 +895,108 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): & (modelled_kwh > floor_kwh) & (measured_kwh >= 0.0) ) - if usable.sum() < 12: + shape_usable = usable.copy() + + # Calibration represents the available PV potential. A battery or inverter + # outage can make an otherwise healthy plant cover only local demand; those + # intervals must not be learned as a permanent model loss. Detect this at day + # level, because individual cloudy hours are much too noisy for a reliable + # availability decision. + local_days = samples_frame.index.tz_convert(self.config.general.timezone).normalize() + fit_start = pd.Timestamp( + end.subtract(days=settings.calibration_days) + .in_timezone(self.config.general.timezone) + .isoformat() + ).normalize() + + if settings.calibration_outage_filter_enabled and usable.any(): + samples = pd.DataFrame( + { + "modelled_kwh": modelled_kwh[usable], + "measured_kwh": measured_kwh[usable], + "local_day": local_days[usable], + } + ) + daily = samples.groupby("local_day").agg( + modelled_kwh=("modelled_kwh", "sum"), + measured_kwh=("measured_kwh", "sum"), + usable_intervals=("modelled_kwh", "size"), + ) + daily["ratio"] = daily["measured_kwh"] / daily["modelled_kwh"] + + # Low-yield weather days do not carry enough evidence to call an outage. + # Half an equivalent full-load hour scales naturally with plant size. + minimum_day_kwh = max(0.5 * self._installed_ac_capacity_w() / 1000.0, 1.0) + reference_candidates = daily[ + (daily["modelled_kwh"] >= minimum_day_kwh) & np.isfinite(daily["ratio"]) + ] + + outage_days = pd.DatetimeIndex([]) + reference_ratio = float("nan") + if len(reference_candidates) >= settings.calibration_min_healthy_days: + # The upper quartile is a robust estimate of the available plant level: + # outages and curtailment only pull the ratio down, while a few weather + # outliers cannot dominate it as a maximum would. + reference_ratio = float(reference_candidates["ratio"].quantile(0.75)) + outage_limit = reference_ratio * settings.calibration_outage_threshold + outage_days = pd.DatetimeIndex( + reference_candidates.index[reference_candidates["ratio"] < outage_limit] + ) + + # A cloudy day inside a known low-production block can accidentally + # resemble a healthy ratio because both numerator and denominator are + # small. Bridge a single-day gap between two detected outage days so a + # continuous plant event is not partly admitted into the fit. + ordered_days = pd.DatetimeIndex(daily.index).sort_values() + bridged_days = [] + for index in range(1, len(ordered_days) - 1): + previous_day = ordered_days[index - 1] + day = ordered_days[index] + next_day = ordered_days[index + 1] + if ( + previous_day in outage_days + and next_day in outage_days + and (day.date() - previous_day.date()).days == 1 + and (next_day.date() - day.date()).days == 1 + ): + bridged_days.append(day) + if bridged_days: + outage_days = outage_days.union(pd.DatetimeIndex(bridged_days)).sort_values() + + healthy_days = pd.DatetimeIndex(daily.index).difference(outage_days).sort_values() + recent_healthy_days = healthy_days[healthy_days >= fit_start] + if len(recent_healthy_days) < settings.calibration_min_healthy_days: + fit_days = healthy_days[-settings.calibration_min_healthy_days :] + else: + fit_days = recent_healthy_days + + # The recent healthy window tracks the current energy level. The stable + # intraday signature uses the full healthy reference window, avoiding noisy + # shape factors learned from only a handful of days. + shape_usable &= np.asarray(local_days.isin(healthy_days), dtype=bool) + usable &= np.asarray(local_days.isin(fit_days), dtype=bool) + if len(outage_days) > 0: + day_list = ", ".join(day.strftime("%Y-%m-%d") for day in outage_days) + logger.info( + "PVForecastAkkudoktorLocal calibration: excluded probable outage/" + f"curtailment days [{day_list}] (healthy reference " + f"{reference_ratio:.3f})." + ) + else: + usable &= np.asarray(local_days >= fit_start, dtype=bool) + shape_usable = usable.copy() + + minimum_samples = math.ceil(12 / interval_hours) + if usable.sum() < minimum_samples: logger.info( - f"PVForecastAkkudoktorLocal calibration: only {int(usable.sum())} usable hours - skipping." + "PVForecastAkkudoktorLocal calibration: only " + f"{int(usable.sum())} usable {interval_minutes}-minute intervals - skipping." ) return None + shape_modelled_kwh = modelled_kwh[shape_usable] + shape_measured_kwh = measured_kwh[shape_usable] + shape_azimuth = azimuth[shape_usable] modelled_kwh = modelled_kwh[usable] measured_kwh = measured_kwh[usable] azimuth = azimuth[usable] @@ -712,54 +1012,83 @@ class PVForecastAkkudoktorLocal(PVForecastProvider): ) ) - bin_degrees = settings.calibration_azimuth_bin_degrees - if bin_degrees <= 0: + factors = self._fit_azimuth_factors( + shape_modelled_kwh, + shape_measured_kwh, + shape_azimuth, + global_factor, + ) + if len(factors) == 1: 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, + f"from {usable.sum()} {interval_minutes}-minute intervals." ) + return global_factor, factors # 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()) + fitted_scale = self._interpolate_azimuth_factors(azimuth, factors) + after = float(np.abs(modelled_kwh * fitted_scale - 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" + f"PVForecastAkkudoktorLocal calibration over {usable.sum()} intervals: global factor " + f"{global_factor:.3f}, {len(factors)} azimuth bins at {interval_minutes}-minute " + "measurement resolution, " + f"MAE {before:.3f} -> {after:.3f} kWh per interval" ) 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.""" + def _interpolate_azimuth_factors(azimuth: np.ndarray, factors: np.ndarray) -> np.ndarray: + """Interpolate fitted bin-centre factors without a discontinuity at north.""" bin_count = len(factors) + if bin_count == 1: + return np.full(len(azimuth), factors[0], dtype=float) + + bin_width = 360.0 / bin_count + centers = (np.arange(bin_count, dtype=float) + 0.5) * bin_width + interpolation_azimuths = np.concatenate( + ([centers[-1] - 360.0], centers, [centers[0] + 360.0]) + ) + interpolation_factors = np.concatenate(([factors[-1]], factors, [factors[0]])) + return np.interp( + np.nan_to_num(azimuth % 360.0), + interpolation_azimuths, + interpolation_factors, + ) + + @staticmethod + def _apply_calibration( + frame: pd.DataFrame, + factors: np.ndarray, + ac_cap_w: float, + global_factor: Optional[float] = None, + timezone: Optional[str] = None, + ) -> pd.DataFrame: + """Scale power with a smooth, circular interpolation of azimuth 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] + scale = PVForecastAkkudoktorLocal._interpolate_azimuth_factors(azimuth, factors) + + # Preserve the independently fitted kWh correction for every forecast day. + # Azimuth factors may redistribute energy within a day, but cannot change its + # calibrated total (apart from the physical inverter cap applied below). + if global_factor is not None and len(frame) > 0: + ac_power = frame["ac_power"].to_numpy(dtype=float) + if isinstance(frame.index, pd.DatetimeIndex): + day_index = frame.index + if timezone is not None and day_index.tz is not None: + day_index = day_index.tz_convert(timezone) + groups = pd.Series(np.arange(len(frame)), index=frame.index).groupby( + day_index.normalize() + ) + group_indices = (group.to_numpy() for _, group in groups) + else: + group_indices = (np.arange(len(frame)),) + for indices in group_indices: + model_energy = float(ac_power[indices].sum()) + shaped_energy = float(np.dot(ac_power[indices], scale[indices])) + if model_energy > 0.0 and shaped_energy > 0.0: + scale[indices] *= global_factor * model_energy / shaped_energy + frame = frame.copy() frame["dc_power"] = frame["dc_power"] * scale frame["ac_power"] = frame["ac_power"] * scale diff --git a/tests/test_measurement.py b/tests/test_measurement.py index 4131f3b4..63f6089b 100644 --- a/tests/test_measurement.py +++ b/tests/test_measurement.py @@ -1,3 +1,5 @@ +import json + import numpy as np import pytest from pendulum import datetime, duration @@ -430,3 +432,30 @@ class TestMeasurement: result = measurement_eos.load_total_kwh(start_datetime=start_datetime, end_datetime=end_datetime, interval=interval) expected = np.array([100]) # Only one complete interval covered np.testing.assert_array_equal(result, expected) + + def test_load_restores_file_records_into_singleton(self, measurement_eos, config_eos, tmp_path): + """File loading must not lose records when validating the Measurement singleton.""" + # ConfigEOS and Measurement are singletons that outlive this module, so + # every global this test touches has to be put back; otherwise later + # modules read their measurements from a deleted tmp_path. + previous_folder = config_eos.general.data_folder_path + previous_keys = config_eos.measurement.load_emr_keys + previous_records = measurement_eos.records + config_eos.general.data_folder_path = tmp_path + config_eos.measurement.load_emr_keys = ["load0_mr"] + record = MeasurementDataRecord(date_time=to_datetime("2026-08-01T12:00:00Z")) + record["load0_mr"] = 123.5 + payload = {"records": [record.model_dump(mode="json")]} + (tmp_path / "measurement.json").write_text( + json.dumps(payload), encoding="utf-8", newline="\n" + ) + + try: + measurement_eos.records = [] + assert measurement_eos.load() is True + assert len(measurement_eos.records) == 1 + assert measurement_eos.records[0]["load0_mr"] == pytest.approx(123.5) + finally: + measurement_eos.records = previous_records + config_eos.measurement.load_emr_keys = previous_keys + config_eos.general.data_folder_path = previous_folder diff --git a/tests/test_pvforecastakkudoktorlocal.py b/tests/test_pvforecastakkudoktorlocal.py index 151c7c69..6af7bf48 100644 --- a/tests/test_pvforecastakkudoktorlocal.py +++ b/tests/test_pvforecastakkudoktorlocal.py @@ -126,7 +126,9 @@ def test_records_are_shifted_to_interval_start(pvforecast_instance): shifted = pvforecast_instance._forecast_frame(synthetic_openmeteo()) pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( - PVForecastAkkudoktorLocalCommonSettings(resolution_minutes=15, shift_to_interval_start=False) + PVForecastAkkudoktorLocalCommonSettings( + resolution_minutes=15, shift_to_interval_start=False + ) ) raw = pvforecast_instance._forecast_frame(synthetic_openmeteo()) @@ -173,7 +175,9 @@ def test_horizon_shading_reduces_yield(pvforecast_instance): def test_update_data_writes_records(pvforecast_instance): - with patch.object(PVForecastAkkudoktorLocal, "_request_forecast", return_value=synthetic_openmeteo()): + with patch.object( + PVForecastAkkudoktorLocal, "_request_forecast", return_value=synthetic_openmeteo() + ): pvforecast_instance._update_data(force_update=True) assert len(pvforecast_instance.records) > 0 @@ -210,6 +214,59 @@ def _feed_measurements( ) +def _feed_measurements_with_recent_outage( + instance: PVForecastAkkudoktorLocal, + frame: pd.DataFrame, + healthy_bias: float, + outage_bias: float, + outage_days: int, + key: str, +) -> None: + """Write a healthy meter history followed by demand-limited PV production.""" + instance.config.measurement.pv_production_emr_keys = [key] + measurement = get_measurement() + + hourly = frame["ac_power"].resample("1h").mean() + hourly = hourly.loc[hourly.index < START] + outage_start = START.subtract(days=outage_days) + healthy_looking_gap = outage_start.add(days=2).date() + + cumulative = 0.0 + for timestamp, power_w in hourly.items(): + measurement.update_value( + pendulum.instance(timestamp.to_pydatetime()), key, round(cumulative, 6) + ) + in_outage = timestamp >= outage_start and timestamp.date() != healthy_looking_gap + bias = outage_bias if in_outage else healthy_bias + cumulative += float(power_w) * bias / 1000.0 + measurement.update_value( + pendulum.instance(hourly.index[-1].to_pydatetime()).add(hours=1), + key, + round(cumulative, 6), + ) + + +def _feed_native_quarter_hour_measurements( + instance: PVForecastAkkudoktorLocal, frame: pd.DataFrame, bias: float, key: str +) -> None: + """Write cumulative PV readings at the provider's native 15-minute cadence.""" + instance.config.measurement.pv_production_emr_keys = [key] + measurement = get_measurement() + slots = frame.loc[frame.index < START, "ac_power"] + + cumulative = 0.0 + for timestamp, power_w in slots.items(): + measurement.update_value( + pendulum.instance(timestamp.to_pydatetime()), key, round(cumulative, 6) + ) + cumulative += float(power_w) * 0.25 * bias / 1000.0 + measurement.update_value( + pendulum.instance(slots.index[-1].to_pydatetime()).add(minutes=15), + 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 @@ -246,6 +303,84 @@ def test_calibration_recovers_a_systematic_bias(pvforecast_instance): assert corrected["ac_power"].sum() == pytest.approx(frame["ac_power"].sum() * global_factor) +def test_calibration_excludes_recent_demand_limited_outage(pvforecast_instance, caplog): + """A battery outage must not teach demand-limited PV as available generation.""" + pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( + PVForecastAkkudoktorLocalCommonSettings( + calibration_enabled=True, + calibration_days=5, + calibration_reference_days=14, + calibration_azimuth_bin_degrees=0, + calibration_min_factor=0.2, + ) + ) + + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + _feed_measurements_with_recent_outage( + pvforecast_instance, + frame, + healthy_bias=0.8, + outage_bias=0.2, + outage_days=5, + key="pv_demand_limited_emr", + ) + + with caplog.at_level("INFO"): + 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]) + assert "excluded probable outage/curtailment days" in caplog.text + # A single statistically healthy-looking day inside the outage is bridged. + assert "2025-06-12" in caplog.text + + # Filtering only selects training data. Applying a global factor preserves the + # native quarter-hour shape instead of replacing it with hourly bucket values. + corrected = pvforecast_instance._apply_calibration(frame, factors, 10000.0) + producing = frame["ac_power"] > 0.0 + assert corrected.index.to_series().diff().dropna().unique().tolist() == [ + pd.Timedelta(minutes=15) + ] + assert ( + corrected.loc[producing, "ac_power"] / frame.loc[producing, "ac_power"] + ).to_numpy() == pytest.approx(np.full(producing.sum(), global_factor)) + + +def test_calibration_uses_real_quarter_hour_measurements(pvforecast_instance, caplog): + """Native meter slots permit shape calibration without inventing intrahour data.""" + pvforecast_instance.config.pvforecast.provider_settings.PVForecastAkkudoktorLocal = ( + PVForecastAkkudoktorLocalCommonSettings( + calibration_enabled=True, + calibration_days=14, + calibration_reference_days=14, + calibration_azimuth_bin_degrees=0, + ) + ) + frame = pvforecast_instance._forecast_frame(synthetic_openmeteo(), calibrate=False) + _feed_native_quarter_hour_measurements( + pvforecast_instance, frame, bias=0.8, key="pv_quarter_hour_emr" + ) + + assert ( + pvforecast_instance._calibration_interval_minutes( + START.subtract(days=14), START + ) + == 15 + ) + with caplog.at_level("INFO"): + calibration = pvforecast_instance._fit_calibration(frame) + + assert calibration is not None + global_factor, _ = calibration + assert global_factor == pytest.approx(0.8, abs=0.03) + assert ( + "15-minute measurement resolution" in caplog.text + or "15-minute intervals" in caplog.text + ) + + 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 = ( @@ -273,6 +408,73 @@ def test_calibration_respects_the_inverter_cap(pvforecast_instance): assert corrected["ac_power"].max() <= 10000.0 + 1e-6 +def test_azimuth_calibration_is_interpolated_smoothly(): + """Azimuth correction must not introduce steps into the quarter-hour plan.""" + frame = pd.DataFrame( + { + "solar_azimuth": [44.9, 45.0, 45.1, 134.9, 135.0, 135.1], + "dc_power": [1000.0] * 6, + "ac_power": [1000.0] * 6, + } + ) + corrected = PVForecastAkkudoktorLocal._apply_calibration( + frame, np.array([0.5, 1.0, 1.5, 1.0]), ac_cap_w=2000.0 + ) + + assert corrected.loc[1, "ac_power"] == pytest.approx(500.0) + assert corrected.loc[4, "ac_power"] == pytest.approx(1000.0) + assert abs(corrected.loc[2, "ac_power"] - corrected.loc[0, "ac_power"]) < 2.0 + assert abs(corrected.loc[5, "ac_power"] - corrected.loc[3, "ac_power"]) < 2.0 + + +def test_azimuth_shape_preserves_each_days_global_energy(): + """The EMS gets a changed shape without a changed daily energy budget.""" + index = pd.date_range("2025-06-01", periods=8, freq="12h", tz="UTC") + frame = pd.DataFrame( + { + "solar_azimuth": [45.0, 225.0, 45.0, 225.0, 45.0, 225.0, 45.0, 225.0], + "dc_power": [100.0, 300.0, 200.0, 200.0, 300.0, 100.0, 150.0, 250.0], + "ac_power": [100.0, 300.0, 200.0, 200.0, 300.0, 100.0, 150.0, 250.0], + }, + index=index, + ) + global_factor = 0.8 + corrected = PVForecastAkkudoktorLocal._apply_calibration( + frame, + np.array([0.5, 1.0, 1.5, 1.0]), + ac_cap_w=10_000.0, + global_factor=global_factor, + timezone="UTC", + ) + + raw_daily = frame["ac_power"].resample("1D").sum() + corrected_daily = corrected["ac_power"].resample("1D").sum() + assert corrected_daily.to_numpy() == pytest.approx( + raw_daily.to_numpy() * global_factor + ) + assert np.std(corrected["ac_power"] / frame["ac_power"]) > 0.01 + + +def test_azimuth_shape_fit_preserves_global_energy(pvforecast_instance): + """Intraday correction must not undo the independently fitted daily kWh.""" + centers = np.arange(22.5, 360.0, 45.0) + azimuth = np.repeat(centers, 20) + modelled_kwh = np.ones(len(azimuth)) + expected_shape = np.repeat([0.8, 0.9, 1.0, 1.1, 1.2, 1.1, 1.0, 0.9], 20) + global_factor = 0.8 + measured_kwh = modelled_kwh * global_factor * expected_shape + + factors = pvforecast_instance._fit_azimuth_factors( + modelled_kwh, measured_kwh, azimuth, global_factor + ) + fitted_scale = pvforecast_instance._interpolate_azimuth_factors(azimuth, factors) + + assert np.std(factors) > 0.01 + assert np.dot(modelled_kwh, fitted_scale) == pytest.approx( + global_factor * modelled_kwh.sum(), rel=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 = (