feat(scada): add calibrated monitoring and incident recovery

scripts/benchmark/calibration.py

@@ -0,0 +1,226 @@
+"""File-isolated hydraulic calibration; database connections are read-only."""
+from __future__ import annotations
+
+import ctypes
+import hashlib
+import json
+import math
+import os
+import re
+import subprocess
+from datetime import datetime, timedelta
+from pathlib import Path
+from uuid import NAMESPACE_URL, uuid5
+from zoneinfo import ZoneInfo
+
+from scripts.benchmark.simulation_design import PipeCalibration, SelectionRule
+
+
+def digest(path: Path) -> str:
+    return hashlib.sha256(path.read_bytes()).hexdigest()
+
+
+def save_json(path: Path, value: dict) -> None:
+    temporary = path.with_suffix(".pending")
+    with temporary.open("x", encoding="utf-8") as output:
+        os.chmod(temporary, 0o600)
+        json.dump(value, output, ensure_ascii=False, indent=2, allow_nan=False, default=str)
+        output.write("\n")
+    temporary.replace(path)
+
+
+def select_calibrated_area(
+    pipe: PipeCalibration, passed_levels: list[str], case_results: list[dict], retain_weak: bool,
+) -> dict:
+    selected = passed_levels[-1] if passed_levels else None
+    soft_blockers = {"insufficient_consecutive_pressure_and_corroboration", "event_rule_did_not_trigger"}
+    physically_valid = bool(case_results) and all(
+        not set(case["blockers"]) - soft_blockers for case in case_results
+    )
+    status = "calibrated_candidate" if selected else (
+        "weak_signal_candidate" if retain_weak and physically_valid else "unsuitable"
+    )
+    return {"status": status, "severity": selected,
+            "area_cm2": pipe.candidate_orifice_areas_cm2[selected] if selected else None}
+
+
+def configure_inp(
+    source: str,
+    start: datetime,
+    duration: int,
+    *,
+    emitter: tuple[str, float] | None = None,
+    pda_options: dict[str, str] | None = None,
+) -> str:
+    """Apply the burst service's options to a file, including the zero-leak control."""
+    if start.tzinfo is None:
+        raise ValueError("simulation start must include a timezone")
+    start = start.astimezone(ZoneInfo("Asia/Shanghai"))
+    updates = {
+        "TIMES": {"DURATION": f"{duration // 3600}:{duration % 3600 // 60:02d}:00",
+                  "PATTERN START": start.strftime("%H:%M:%S"),
+                  "REPORT START": "0:00:00", "STATISTIC": "NONE"},
+        "OPTIONS": {**(pda_options or {}), "DEMAND MODEL": "PDA", "REQUIRED PRESSURE": "10.0000"},
+        "REPORT": {"NODES": "NONE", "LINKS": "NONE"},
+    }
+    section = None
+    output = []
+    pending = {}
+    for line in source.splitlines():
+        match = re.fullmatch(r"\s*\[([^]]+)\]\s*", line)
+        if match:
+            output.extend(f"{key} {value}" for key, value in pending.items())
+            section = match[1].upper()
+            pending = dict(updates.get(section, {}))
+            output.append(line)
+            if section == "EMITTERS" and emitter:
+                node, area = emitter
+                coefficient = 0.65 * area / 10000 * math.sqrt(19.6) * 1000
+                output.append(f"{node} {coefficient:.12g}")
+            continue
+        tokens = line.split(";", 1)[0].upper().split()
+        if any(tokens[:len(key.split())] == key.split() for key in updates.get(section, {})):
+            continue
+        output.append(line)
+    output.extend(f"{key} {value}" for key, value in pending.items())
+    if "[EMITTERS]" not in source.upper() and emitter:
+        raise ValueError("exported model has no emitter section")
+    return "\n".join(output) + "\n"
+
+
+def run_file(inp: Path, node_ids: list[str], link_ids: list[str], *, timeout: int) -> dict:
+    from app.infra.epanet.epanet import Output
+
+    executable = Path(__file__).resolve().parents[2] / "app/infra/epanet/linux/runepanet"
+    env = os.environ.copy()
+    env["LD_LIBRARY_PATH"] = f"{executable.parent}:{env.get('LD_LIBRARY_PATH', '')}"
+    report_path, binary = inp.with_suffix(".rpt"), inp.with_suffix(".bin")
+    result = subprocess.run(
+        [str(executable), str(inp), str(report_path), str(binary)],
+        env=env, capture_output=True, text=True, timeout=timeout, check=False,
+    )
+    report = report_path.read_text(encoding="utf-8", errors="replace") if report_path.exists() else ""
+    warnings = [line.strip() for line in report.splitlines()
+                if re.search(r"^\s*(?:warning\b|error\s+\d+\b)|system\s+(?:is\s+)?unbalanced|cannot solve hydraulic",
+                             line, re.IGNORECASE)]
+    summary = {"exit_code": result.returncode, "warnings": warnings,
+               "inp_sha256": digest(inp), "report_sha256": digest(report_path) if report_path.exists() else None}
+    if result.returncode or not binary.exists():
+        return {**summary, "status": "failed"}
+    reader = Output(str(binary))
+    names = reader.element_name()
+    summary.update(status="completed", units=reader.units(), times=reader.times(), nodes={}, links={})
+    for kind, selected, attributes, function in (
+        ("nodes", node_ids, ["demand", "head", "pressure", "quality"], reader._lib.ENR_getNodeResult),
+        ("links", link_ids, ["flow", "velocity", "headloss", "quality", "status", "setting", "reaction", "friction"],
+         reader._lib.ENR_getLinkResult),
+    ):
+        indices = {name: index + 1 for index, name in enumerate(names[kind])}
+        for name in selected:
+            index = indices.get(name)
+            if index is None:
+                raise ValueError("requested hydraulic output is missing")
+            series = []
+            for period in range(summary["times"]["num_periods"]):
+                values = ctypes.POINTER(ctypes.c_float)()
+                length = ctypes.c_int()
+                reader._check(function(reader._handle, period, index, ctypes.byref(values), ctypes.byref(length)))
+                try:
+                    if length.value != len(attributes):
+                        raise ValueError("unexpected hydraulic result attributes")
+                    row = {key: float(values[i]) for i, key in enumerate(attributes)}
+                    if not all(math.isfinite(value) for value in row.values()):
+                        raise ValueError("non-finite hydraulic result")
+                    series.append(row)
+                finally:
+                    reader._lib.ENR_free(ctypes.byref(values))
+            summary[kind][name] = series
+    return summary
+
+
+def pressure_preview(
+    pipe: PipeCalibration, start: datetime, duration: int, step: int,
+    history: dict, current: dict, rule: SelectionRule,
+) -> dict:
+    """Exercise the actual replay cleaner/detector without database mutations."""
+    from app.domain.schemas.scada_replay import ScadaReplayCreate
+    from app.services.scada_replay_processing import detect_frame, prepare_cleaning
+
+    ids = [binding.device_id for binding in pipe.pressure_bindings]
+    config = ScadaReplayCreate(
+        request_id=uuid5(NAMESPACE_URL, f"calibration:{pipe.pipe_id}:{start.isoformat()}"),
+        start_time=start, end_time=start + timedelta(seconds=duration + step), device_ids=ids,
+        rules=[{"rule_id": "calibration-pressure", "device_ids": ids,
+                "kind": "relative_change", "direction": "decrease", "unit": "m",
+                "threshold": rule.primary_pressure_drop_percent,
+                "window_samples": 4, "minimum_samples": 4,
+                "consecutive_samples": rule.consecutive_samples, "minimum_devices": 1}],
+    )
+    document = {"replay_id": str(config.request_id), "devices": [
+        {"device_id": key, "transmission_frequency": "0:15:00"} for key in ids],
+        "histories": {binding.device_id: [row["pressure"] for row in history[binding.node_id]]
+                      for binding in pipe.pressure_bindings}}
+    raised = []
+    for index in range(duration // step + 1):
+        at = start + timedelta(seconds=index * step)
+        rows = [{"device_id": binding.device_id, "time": at,
+                 "monitored_value": current[binding.node_id][index]["pressure"],
+                 "source_hash": hashlib.sha256(f"{binding.device_id}:{at}".encode()).hexdigest()}
+                for binding in pipe.pressure_bindings]
+        document, changes = detect_frame(document, prepare_cleaning(rows, config), config, at)
+        raised.extend(change for change in changes
+                      if change["change"] == "raised" and change.get("category") != "data_quality")
+    return {"raised": bool(raised), "first_signal_at": raised[0]["occurred_at"] if raised else None,
+            "scope": "pure replay rule preview; not monitoring worker or Agent acceptance"}
+
+
+def assess(
+    pipe: PipeCalibration, normal: dict, burst: dict, rule: SelectionRule,
+    *, duration: int, step: int,
+) -> dict:
+    expected_times = {"report_start": 0, "report_step": step,
+                      "sim_duration": duration, "num_periods": duration // step + 1}
+    blockers = []
+    for label, result in (("normal", normal), ("burst", burst)):
+        if result.get("status") != "completed" or result.get("times") != expected_times:
+            blockers.append(f"{label}_incomplete_or_failed")
+        if result.get("warnings"):
+            blockers.append(f"{label}_solver_warning")
+        if result.get("units") != {"flow": "LPS", "pressure": "MTR", "quality": "NONE"}:
+            blockers.append(f"{label}_unexpected_units")
+    if blockers:
+        return {"passed": False, "blockers": blockers}
+    pressures, flows = {}, {}
+    for binding in pipe.pressure_bindings:
+        before = [row["pressure"] for row in normal["nodes"][binding.node_id]]
+        after = [row["pressure"] for row in burst["nodes"][binding.node_id]]
+        if min(before) <= rule.relative_change_epsilon:
+            blockers.append(f"nonpositive_normal_pressure:{binding.device_id}")
+        if min(after) < 0:
+            blockers.append(f"negative_monitor_pressure:{binding.device_id}")
+        pressures[binding.device_id] = [100 * (a - b) / max(abs(a), rule.relative_change_epsilon)
+                                        for a, b in zip(before, after)]
+    for binding in pipe.flow_bindings:
+        flows[binding.device_id] = [100 * abs(a["flow"] - b["flow"]) /
+                                   max(abs(a["flow"]), rule.relative_change_epsilon)
+                                   for a, b in zip(normal["links"][binding.link_id], burst["links"][binding.link_id])]
+    witnesses = []
+    count = duration // step + 1
+    for primary, drops in pressures.items():
+        for start in range(count - rule.consecutive_samples + 1):
+            window = range(start, start + rule.consecutive_samples)
+            if not all(drops[i] >= rule.primary_pressure_drop_percent for i in window):
+                continue
+            second = [key for key, values in pressures.items() if key != primary and
+                      all(values[i] >= rule.corroborating_pressure_drop_percent for i in window)]
+            flow = [key for key, values in flows.items() if
+                    all(values[i] >= rule.flow_change_percent for i in window) and
+                    all(abs(normal["links"][next(b.link_id for b in pipe.flow_bindings if b.device_id == key)][i]["flow"])
+                        > rule.relative_change_epsilon for i in window)]
+            if second or flow:
+                witnesses.append({"primary_device": primary, "second_pressure_devices": second,
+                                  "flow_devices": flow, "start_offset_seconds": start * step})
+    if not witnesses:
+        blockers.append("insufficient_consecutive_pressure_and_corroboration")
+    return {"passed": not blockers, "blockers": blockers, "witnesses": witnesses,
+            "pressure_drop_percent": pressures, "flow_change_percent": flows}