Feat/dynamic spatial temporal risk

benchmarks/dynamic/level2/SPATIAL_TEMPORAL_RISK.md

@@ -0,0 +1,22 @@
+# Spatial-temporal voltage risk metrics
+
+The current transient-stability KPI averages normalized voltage excursions
+across buses. That system-level scalar can hide a severe excursion at one bus
+or a brief extreme event. This extension reports spatial and temporal
+diagnostics while leaving the existing severity score unchanged.
+
+Under-voltage is normalized as `clip((0.95 - V)/(0.95 - 0.90), 0, 1)`;
+over-voltage as `clip((V - 1.05)/(1.10 - 1.05), 0, 1)`. Metrics include
+trapezoidally integrated mean and worst-bus exposures, the fraction of buses
+that ever violate a limit, an upper-tail mean of sample-level spatial-average
+violations, and the longest contiguous interval where both interval endpoints
+violate a threshold. Nonuniform time sampling is supported. Duplicate or
+decreasing timestamps and non-finite values are rejected.
+
+The upper-tail mean is CVaR-like in spirit, but is not formal CVaR over a
+probability distribution. These descriptive diagnostics do not replace
+protection limits, stability criteria, or the generator-speed instability flag.
+
+## Offline validation
+Run `python -m unittest -v test_spatial_temporal_risk.py` from this directory.
+Tests use synthetic trajectories and need neither PowerFactory nor network access.

benchmarks/dynamic/level2/spatial_temporal_risk.py

@@ -0,0 +1,85 @@
+"""Additional risk metrics for transient-stability voltage trajectories.
+These metrics complement the existing Integral of Violation (IoV) score.
+They quantify spatial concentration and temporal persistence rather than
+averaging every bus and time sample into one scalar.
+"""
+from __future__ import annotations
+from typing import Any
+import numpy as np
+import pandas as pd
+
+def spatial_temporal_risk(voltages: pd.DataFrame, times: pd.Series, *,
+    undervoltage_threshold: float = 0.95, undervoltage_floor: float = 0.90,
+    overvoltage_threshold: float = 1.05, overvoltage_ceiling: float = 1.10,
+    tail_fraction: float = 0.10) -> dict[str, Any]:
+    """Compute exposure, concentration, tail, and persistence diagnostics.
+
+    Rows are time samples, columns are buses, and time units are seconds.
+    Results are descriptive diagnostics, not protection criteria.
+    """
+    if not 0.0 < tail_fraction <= 1.0:
+        raise ValueError("tail_fraction must be in (0, 1]")
+    if not 0.0 <= undervoltage_floor < undervoltage_threshold:
+        raise ValueError("undervoltage_floor must be below undervoltage_threshold")
+    if not 0.0 < overvoltage_threshold < overvoltage_ceiling:
+        raise ValueError("overvoltage_threshold must be below overvoltage_ceiling")
+    if len(voltages) != len(times):
+        raise ValueError("voltages and times must have the same number of samples")
+    if voltages.empty or voltages.shape[1] == 0:
+        return _empty_result(tail_fraction)
+    values = voltages.to_numpy(dtype=float)
+    if not np.isfinite(values).all():
+        raise ValueError("voltages must contain only finite numeric values")
+    t = np.asarray(times, dtype=float)
+    if not np.isfinite(t).all():
+        raise ValueError("times must contain only finite numeric values")
+    if len(t) > 1 and np.any(np.diff(t) <= 0):
+        raise ValueError("times must be strictly increasing")
+    window = float(t[-1] - t[0]) if len(t) > 1 else 0.0
+    under = np.clip((undervoltage_threshold - values) / (undervoltage_threshold - undervoltage_floor), 0.0, 1.0)
+    over = np.clip((values - overvoltage_threshold) / (overvoltage_ceiling - overvoltage_threshold), 0.0, 1.0)
+    if len(t) < 2 or window <= 0:
+        under_iov = np.zeros(values.shape[1]); over_iov = np.zeros(values.shape[1])
+        under_sample = np.zeros(len(t)); over_sample = np.zeros(len(t))
+        under_duration = np.zeros(values.shape[1]); over_duration = np.zeros(values.shape[1])
+    else:
+        # Trapezoidal integration accounts for nonuniform solver sampling.
+        under_iov = np.trapezoid(under, t, axis=0) / window
+        over_iov = np.trapezoid(over, t, axis=0) / window
+        under_sample = under.mean(axis=1); over_sample = over.mean(axis=1)
+        intervals = np.diff(t)
+        under_duration = ((under[:-1] > 0) & (under[1:] > 0)).T @ intervals
+        over_duration = ((over[:-1] > 0) & (over[1:] > 0)).T @ intervals
+    return {
+        "sample_count": int(len(t)), "bus_count": int(values.shape[1]),
+        "window_seconds": _r(window),
+        "mean_iov_undervoltage": _r(float(under_iov.mean())),
+        "mean_iov_overvoltage": _r(float(over_iov.mean())),
+        "worst_bus_iov_undervoltage": _r(float(under_iov.max())),
+        "worst_bus_iov_overvoltage": _r(float(over_iov.max())),
+        "affected_bus_fraction_undervoltage": _r(float((under > 0).any(axis=0).mean())),
+        "affected_bus_fraction_overvoltage": _r(float((over > 0).any(axis=0).mean())),
+        "tail_mean_undervoltage": _r(_tail_mean(under_sample, tail_fraction)),
+        "tail_mean_overvoltage": _r(_tail_mean(over_sample, tail_fraction)),
+        "longest_violation_seconds_undervoltage": _r(float(under_duration.max(initial=0.0))),
+        "longest_violation_seconds_overvoltage": _r(float(over_duration.max(initial=0.0))),
+        "tail_fraction": float(tail_fraction),
+    }
+
+def _tail_mean(values: np.ndarray, fraction: float) -> float:
+    if values.size == 0: return 0.0
+    count = max(1, int(np.ceil(values.size * fraction)))
+    return float(np.sort(values)[-count:].mean())
+
+def _r(value: float) -> float:
+    return round(value, 8)
+
+def _empty_result(fraction: float) -> dict[str, Any]:
+    keys = ("mean_iov_undervoltage", "mean_iov_overvoltage",
+        "worst_bus_iov_undervoltage", "worst_bus_iov_overvoltage",
+        "affected_bus_fraction_undervoltage", "affected_bus_fraction_overvoltage",
+        "tail_mean_undervoltage", "tail_mean_overvoltage",
+        "longest_violation_seconds_undervoltage", "longest_violation_seconds_overvoltage")
+    result = {key: 0.0 for key in keys}
+    result.update(sample_count=0, bus_count=0, window_seconds=0.0, tail_fraction=float(fraction))
+    return result

benchmarks/dynamic/level2/test_spatial_temporal_risk.py

@@ -0,0 +1,51 @@
+import unittest
+import numpy as np
+import pandas as pd
+from spatial_temporal_risk import spatial_temporal_risk
+
+class SpatialTemporalRiskTests(unittest.TestCase):
+    def test_localized_and_widespread_disturbance_are_distinguished(self):
+        t = pd.Series([0.0, 1.0, 2.0])
+        local = spatial_temporal_risk(pd.DataFrame({"A": [0.90]*3, "B": [1.0]*3}), t)
+        wide = spatial_temporal_risk(pd.DataFrame({"A": [0.90]*3, "B": [0.90]*3}), t)
+        self.assertLess(local["mean_iov_undervoltage"], wide["mean_iov_undervoltage"])
+        self.assertEqual(local["affected_bus_fraction_undervoltage"], 0.5)
+        self.assertEqual(wide["affected_bus_fraction_undervoltage"], 1.0)
+
+    def test_tail_detects_brief_severe_excursion(self):
+        t = pd.Series([0., 1., 2., 3., 4.])
+        v = pd.DataFrame({"A": [1., 1., .90, 1., 1.], "B": [1.]*5})
+        result = spatial_temporal_risk(v, t, tail_fraction=.2)
+        self.assertGreater(result["tail_mean_undervoltage"], result["mean_iov_undervoltage"])
+
+    def test_irregular_sample_times(self):
+        result = spatial_temporal_risk(pd.DataFrame({"A": [.95, .90, .95]}), pd.Series([0., 1., 3.]))
+        self.assertAlmostEqual(result["mean_iov_undervoltage"], .5)
+
+    def test_under_over_metrics_are_independent(self):
+        t = pd.Series([0., 1., 2.])
+        under = spatial_temporal_risk(pd.DataFrame({"A": [.90]*3}), t)
+        over = spatial_temporal_risk(pd.DataFrame({"A": [1.10]*3}), t)
+        self.assertEqual(under["worst_bus_iov_overvoltage"], 0.0)
+        self.assertEqual(over["worst_bus_iov_undervoltage"], 0.0)
+        self.assertEqual(over["worst_bus_iov_overvoltage"], 1.0)
+
+    def test_empty_and_single_sample_inputs(self):
+        empty = spatial_temporal_risk(pd.DataFrame(), pd.Series(dtype=float))
+        one = spatial_temporal_risk(pd.DataFrame({"A": [.90]}), pd.Series([0.]))
+        self.assertEqual(empty["bus_count"], 0)
+        self.assertEqual(one["mean_iov_undervoltage"], 0.0)
+
+    def test_invalid_contract(self):
+        v = pd.DataFrame({"A": [1., 1.]})
+        with self.assertRaisesRegex(ValueError, "same number"):
+            spatial_temporal_risk(v, pd.Series([0.]))
+        with self.assertRaisesRegex(ValueError, "strictly increasing"):
+            spatial_temporal_risk(v, pd.Series([1., 1.]))
+        with self.assertRaisesRegex(ValueError, "finite"):
+            spatial_temporal_risk(pd.DataFrame({"A": [1., np.nan]}), pd.Series([0., 1.]))
+        with self.assertRaisesRegex(ValueError, "tail_fraction"):
+            spatial_temporal_risk(v, pd.Series([0., 1.]), tail_fraction=0.)
+
+if __name__ == "__main__":
+    unittest.main()