feat(benchmark): add tjwater2 burst calibration design
This commit is contained in:
@@ -0,0 +1,185 @@
|
||||
"""Strict contracts for a burst-simulation calibration design."""
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from datetime import datetime
|
||||
from pathlib import Path
|
||||
from typing import Literal
|
||||
from uuid import UUID
|
||||
|
||||
from pydantic import BaseModel, ConfigDict, Field, model_validator
|
||||
|
||||
|
||||
Severity = Literal["light", "medium", "severe"]
|
||||
|
||||
|
||||
class StrictModel(BaseModel):
|
||||
model_config = ConfigDict(extra="forbid", allow_inf_nan=False)
|
||||
|
||||
|
||||
class ModelExpectations(StrictModel):
|
||||
node_count: int = Field(gt=0)
|
||||
pipe_count: int = Field(gt=0)
|
||||
valve_count: int = Field(gt=0)
|
||||
flow_unit: Literal["LPS"]
|
||||
pressure_unit: Literal["METERS"]
|
||||
pressure_measurement_unit: Literal["m"]
|
||||
flow_measurement_unit: Literal["m3/h"]
|
||||
emitter_exponent: float = Field(gt=0)
|
||||
hydraulic_timestep_seconds: int = Field(gt=0)
|
||||
report_timestep_seconds: int = Field(gt=0)
|
||||
|
||||
|
||||
class SeverityLevel(StrictModel):
|
||||
level: Severity
|
||||
cross_section_fraction: float = Field(gt=0, le=0.125)
|
||||
|
||||
|
||||
class SelectionRule(StrictModel):
|
||||
primary_pressure_drop_percent: float = Field(gt=0, le=100)
|
||||
corroborating_pressure_drop_percent: float = Field(gt=0, le=100)
|
||||
flow_change_percent: float = Field(gt=0, le=1000)
|
||||
consecutive_samples: int = Field(ge=1, le=20)
|
||||
expected_sample_count: int = Field(ge=2, le=100)
|
||||
relative_change_epsilon: float = Field(gt=0, le=0.001)
|
||||
normal_control_must_not_trigger: Literal[True]
|
||||
negative_monitor_pressure_allowed: Literal[False]
|
||||
choose_smallest_passing_area: Literal[True]
|
||||
|
||||
|
||||
class BindingValidation(StrictModel):
|
||||
max_pressure_distance_m: float = Field(gt=0, le=10)
|
||||
max_flow_distance_m: float = Field(gt=0, le=10)
|
||||
|
||||
|
||||
class SimulationCase(StrictModel):
|
||||
case_id: str = Field(pattern=r"^B[0-9]{2}$")
|
||||
pipe_id: str = Field(min_length=1)
|
||||
occurred_at: datetime
|
||||
|
||||
@model_validator(mode="after")
|
||||
def aware_time(self):
|
||||
if self.occurred_at.tzinfo is None:
|
||||
raise ValueError("case time must include a timezone")
|
||||
return self
|
||||
|
||||
|
||||
class PressureBinding(StrictModel):
|
||||
device_id: str = Field(min_length=1)
|
||||
node_id: str = Field(min_length=1)
|
||||
|
||||
|
||||
class FlowBinding(StrictModel):
|
||||
device_id: str = Field(min_length=1)
|
||||
link_id: str = Field(min_length=1)
|
||||
|
||||
|
||||
class PipeCalibration(StrictModel):
|
||||
role: Literal["primary", "secondary"]
|
||||
pipe_id: str = Field(min_length=1)
|
||||
diameter_mm: float = Field(gt=0)
|
||||
length_m: float = Field(gt=0)
|
||||
start_node_id: str = Field(min_length=1)
|
||||
emitter_node_id: str = Field(min_length=1)
|
||||
longitude: float = Field(ge=-180, le=180)
|
||||
latitude: float = Field(gt=-90, lt=90)
|
||||
candidate_orifice_areas_cm2: dict[Severity, float]
|
||||
pressure_bindings: list[PressureBinding] = Field(min_length=1, max_length=20)
|
||||
flow_bindings: list[FlowBinding] = Field(min_length=1, max_length=20)
|
||||
boundary_valve_ids: list[str] = Field(min_length=1, max_length=20)
|
||||
output_node_ids: list[str] = Field(min_length=1, max_length=20)
|
||||
output_link_ids: list[str] = Field(default_factory=list, max_length=20)
|
||||
|
||||
@model_validator(mode="after")
|
||||
def complete_outputs(self):
|
||||
if set(self.candidate_orifice_areas_cm2) != {"light", "medium", "severe"}:
|
||||
raise ValueError("all three severity areas are required")
|
||||
pressure_device_ids = [value.device_id for value in self.pressure_bindings]
|
||||
pressure_node_ids = [value.node_id for value in self.pressure_bindings]
|
||||
flow_device_ids = [value.device_id for value in self.flow_bindings]
|
||||
flow_link_ids = [value.link_id for value in self.flow_bindings]
|
||||
for values, label in (
|
||||
(pressure_device_ids, "pressure devices"),
|
||||
(pressure_node_ids, "pressure nodes"),
|
||||
(flow_device_ids, "flow devices"),
|
||||
(flow_link_ids, "flow links"),
|
||||
(self.boundary_valve_ids, "boundary valves"),
|
||||
(self.output_node_ids, "output nodes"),
|
||||
(self.output_link_ids, "output links"),
|
||||
):
|
||||
if len(values) != len(set(values)):
|
||||
raise ValueError(f"{label} must be distinct")
|
||||
if not {self.start_node_id, self.emitter_node_id, *pressure_node_ids}.issubset(
|
||||
self.output_node_ids
|
||||
):
|
||||
raise ValueError("pipe endpoints and pressure nodes must be returned")
|
||||
if not {self.pipe_id, *flow_link_ids, *self.boundary_valve_ids}.issubset(
|
||||
self.output_link_ids
|
||||
):
|
||||
raise ValueError("target, flow and boundary links must be returned")
|
||||
return self
|
||||
|
||||
|
||||
class BurstSimulationDesign(StrictModel):
|
||||
version: Literal["burst-simulation-design@1"]
|
||||
suite_id: Literal["tjwater2-v2-burst-calibration"]
|
||||
project_code: Literal["tjwater2_v2"]
|
||||
project_id: UUID
|
||||
model_expectations: ModelExpectations
|
||||
duration_seconds: int = Field(ge=60, le=86400)
|
||||
event_times: list[datetime] = Field(min_length=1, max_length=8)
|
||||
severity_levels: list[SeverityLevel] = Field(min_length=3, max_length=3)
|
||||
calibration_order: list[Severity] = Field(min_length=3, max_length=3)
|
||||
selection_rule: SelectionRule
|
||||
binding_validation: BindingValidation
|
||||
cases: list[SimulationCase] = Field(min_length=1, max_length=20)
|
||||
pipes: list[PipeCalibration] = Field(min_length=1, max_length=4)
|
||||
|
||||
@model_validator(mode="after")
|
||||
def coherent_design(self):
|
||||
if any(value.tzinfo is None for value in self.event_times):
|
||||
raise ValueError("event times must include a timezone")
|
||||
if len(set(self.event_times)) != len(self.event_times):
|
||||
raise ValueError("event times must be distinct")
|
||||
levels = {entry.level: entry.cross_section_fraction for entry in self.severity_levels}
|
||||
if set(levels) != {"light", "medium", "severe"}:
|
||||
raise ValueError("light, medium and severe levels are required")
|
||||
if not levels["light"] < levels["medium"] < levels["severe"]:
|
||||
raise ValueError("severity fractions must increase")
|
||||
expected_levels = {"light": 1 / 32, "medium": 1 / 16, "severe": 1 / 8}
|
||||
if any(
|
||||
not math.isclose(levels[level], fraction, rel_tol=0, abs_tol=1e-12)
|
||||
for level, fraction in expected_levels.items()
|
||||
):
|
||||
raise ValueError("severity fractions must remain 1/32, 1/16 and 1/8")
|
||||
if self.calibration_order != ["severe", "medium", "light"]:
|
||||
raise ValueError("calibration order must be severe, medium, light")
|
||||
report_step = self.model_expectations.report_timestep_seconds
|
||||
if self.duration_seconds % report_step:
|
||||
raise ValueError("duration must contain whole report timesteps")
|
||||
expected_samples = self.duration_seconds // report_step + 1
|
||||
if self.selection_rule.expected_sample_count != expected_samples:
|
||||
raise ValueError("expected sample count does not match duration and report step")
|
||||
pipes = {pipe.pipe_id: pipe for pipe in self.pipes}
|
||||
if len(pipes) != len(self.pipes):
|
||||
raise ValueError("pipe IDs must be distinct")
|
||||
if len({case.case_id for case in self.cases}) != len(self.cases):
|
||||
raise ValueError("case IDs must be distinct")
|
||||
expected_cases = {(pipe_id, at) for pipe_id in pipes for at in self.event_times}
|
||||
actual_cases = {(case.pipe_id, case.occurred_at) for case in self.cases}
|
||||
if actual_cases != expected_cases:
|
||||
raise ValueError("cases must cover every pipe and event time exactly once")
|
||||
for pipe in self.pipes:
|
||||
cross_section_cm2 = math.pi * (pipe.diameter_mm / 10) ** 2 / 4
|
||||
for level, fraction in levels.items():
|
||||
expected = cross_section_cm2 * fraction
|
||||
actual = pipe.candidate_orifice_areas_cm2[level]
|
||||
if not math.isclose(actual, expected, rel_tol=0, abs_tol=0.001):
|
||||
raise ValueError(
|
||||
f"{pipe.pipe_id} {level} area does not match its diameter fraction"
|
||||
)
|
||||
return self
|
||||
|
||||
|
||||
def load_simulation_design(path: Path) -> BurstSimulationDesign:
|
||||
return BurstSimulationDesign.model_validate_json(path.read_text(encoding="utf-8"))
|
||||
@@ -0,0 +1,288 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Read-only validation of the tjwater2_v2 burst-simulation design."""
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import asyncio
|
||||
import sys
|
||||
import traceback
|
||||
from datetime import timedelta
|
||||
from pathlib import Path
|
||||
|
||||
from psycopg import connect
|
||||
from psycopg.conninfo import conninfo_to_dict, make_conninfo
|
||||
from psycopg.rows import dict_row
|
||||
|
||||
|
||||
def readonly(dsn: str):
|
||||
return connect(
|
||||
make_conninfo(
|
||||
dsn,
|
||||
options="-c default_transaction_read_only=on -c statement_timeout=20000",
|
||||
application_name="tjwater-burst-design-readonly",
|
||||
connect_timeout="5",
|
||||
),
|
||||
row_factory=dict_row,
|
||||
)
|
||||
|
||||
|
||||
async def resolve_project(design):
|
||||
from sqlalchemy import text
|
||||
|
||||
from app.infra.db.metadb.database import SessionLocal, engine
|
||||
from app.infra.db.metadb.repositories.metadata_repository import MetadataRepository
|
||||
|
||||
try:
|
||||
async with SessionLocal() as session:
|
||||
await session.execute(
|
||||
text("SET TRANSACTION ISOLATION LEVEL REPEATABLE READ, READ ONLY")
|
||||
)
|
||||
repo = MetadataRepository(session)
|
||||
project = await repo.get_project_by_code(design.project_code)
|
||||
if (
|
||||
project is None
|
||||
or project.id != design.project_id
|
||||
or project.status != "active"
|
||||
):
|
||||
raise ValueError("project identity or active status does not match")
|
||||
biz = await repo.get_project_db_routing(project.id, "biz_data")
|
||||
iot = await repo.get_project_db_routing(project.id, "iot_data")
|
||||
if biz is None or iot is None:
|
||||
raise ValueError("project database routing is incomplete")
|
||||
await session.rollback()
|
||||
return biz.dsn, iot.dsn
|
||||
finally:
|
||||
await engine.dispose()
|
||||
|
||||
|
||||
def validate_database(connection, design, *, template: bool) -> None:
|
||||
connection.execute("SET TRANSACTION ISOLATION LEVEL REPEATABLE READ, READ ONLY")
|
||||
if connection.execute("SHOW transaction_read_only").fetchone()["transaction_read_only"] != "on":
|
||||
raise ValueError("database transaction is not read-only")
|
||||
expected = design.model_expectations
|
||||
counts = connection.execute(
|
||||
"SELECT (SELECT count(*) FROM network.nodes) AS nodes, "
|
||||
"(SELECT count(*) FROM network.pipes) AS pipes, "
|
||||
"(SELECT count(*) FROM network.links WHERE link_type='valve') AS valves"
|
||||
).fetchone()
|
||||
if (counts["nodes"], counts["pipes"], counts["valves"]) != (
|
||||
expected.node_count,
|
||||
expected.pipe_count,
|
||||
expected.valve_count,
|
||||
):
|
||||
raise ValueError("model counts do not match the frozen design")
|
||||
|
||||
settings = {
|
||||
row["key"]: row["value"]
|
||||
for row in connection.execute(
|
||||
"SELECT key,value FROM network.simulation_settings "
|
||||
"WHERE engine_version='legacy' AND key=ANY(%s)",
|
||||
(["UNITS", "PRESSURE", "EMITTER EXPONENT"],),
|
||||
)
|
||||
}
|
||||
times = {
|
||||
row["key"]: row["value"]
|
||||
for row in connection.execute(
|
||||
"SELECT key,value FROM network.time_settings "
|
||||
"WHERE key=ANY(%s)",
|
||||
(["HYDRAULIC TIMESTEP", "REPORT TIMESTEP"],),
|
||||
)
|
||||
}
|
||||
if settings != {
|
||||
"UNITS": expected.flow_unit,
|
||||
"PRESSURE": expected.pressure_unit,
|
||||
"EMITTER EXPONENT": f"{expected.emitter_exponent:.4f}",
|
||||
}:
|
||||
raise ValueError("simulation settings do not match the frozen design")
|
||||
if times != {"HYDRAULIC TIMESTEP": "0:15:00", "REPORT TIMESTEP": "0:15:00"}:
|
||||
raise ValueError("simulation timesteps do not match the frozen design")
|
||||
|
||||
for pipe in design.pipes:
|
||||
row = connection.execute(
|
||||
"SELECT l.start_node_id,l.end_node_id,p.length,p.diameter,p.status "
|
||||
"FROM network.links l JOIN network.pipes p ON p.link_id=l.id WHERE l.id=%s",
|
||||
(pipe.pipe_id,),
|
||||
).fetchone()
|
||||
if row is None:
|
||||
raise ValueError(f"pipe {pipe.pipe_id} is missing")
|
||||
if (
|
||||
str(row["start_node_id"]) != pipe.start_node_id
|
||||
or str(row["end_node_id"]) != pipe.emitter_node_id
|
||||
or abs(float(row["length"]) - pipe.length_m) > 0.001
|
||||
or abs(float(row["diameter"]) - pipe.diameter_mm) > 0.001
|
||||
or str(row["status"]).upper() != "OPEN"
|
||||
):
|
||||
raise ValueError(f"pipe {pipe.pipe_id} properties changed")
|
||||
node_rows = connection.execute(
|
||||
"SELECT id,node_type FROM network.nodes WHERE id=ANY(%s)",
|
||||
(pipe.output_node_ids,),
|
||||
).fetchall()
|
||||
if {str(row["id"]) for row in node_rows} != set(pipe.output_node_ids):
|
||||
raise ValueError(f"pipe {pipe.pipe_id} output nodes are incomplete")
|
||||
node_types = {str(row["id"]): str(row["node_type"]) for row in node_rows}
|
||||
if node_types.get(pipe.emitter_node_id) != "junction":
|
||||
raise ValueError(f"pipe {pipe.pipe_id} emitter endpoint is not a junction")
|
||||
link_rows = connection.execute(
|
||||
"SELECT id,link_type FROM network.links WHERE id=ANY(%s)",
|
||||
(pipe.output_link_ids,),
|
||||
).fetchall()
|
||||
if {str(row["id"]) for row in link_rows} != set(pipe.output_link_ids):
|
||||
raise ValueError(f"pipe {pipe.pipe_id} output links are incomplete")
|
||||
link_types = {str(row["id"]): str(row["link_type"]) for row in link_rows}
|
||||
if any(link_types.get(value) != "valve" for value in pipe.boundary_valve_ids):
|
||||
raise ValueError(f"pipe {pipe.pipe_id} boundary valve changed")
|
||||
emitter = connection.execute(
|
||||
"SELECT coefficient FROM network.emitters WHERE junction_id=%s",
|
||||
(pipe.emitter_node_id,),
|
||||
).fetchone()
|
||||
if emitter is not None:
|
||||
raise ValueError(f"pipe {pipe.pipe_id} endpoint already has an emitter")
|
||||
|
||||
if not template:
|
||||
pressure_bindings = {
|
||||
(
|
||||
str(row["device_id"]),
|
||||
str(row["node_id"]),
|
||||
str(row["measurement_unit"]),
|
||||
str(row["transmission_frequency"]),
|
||||
)
|
||||
for row in connection.execute(
|
||||
"SELECT device_id,node_id,measurement_unit,transmission_frequency "
|
||||
"FROM asset.scada_devices "
|
||||
"WHERE device_type='pressure' AND device_id=ANY(%s)",
|
||||
([value.device_id for value in pipe.pressure_bindings],),
|
||||
)
|
||||
}
|
||||
flow_bindings = {
|
||||
(
|
||||
str(row["device_id"]),
|
||||
str(row["link_id"]),
|
||||
str(row["measurement_unit"]),
|
||||
str(row["transmission_frequency"]),
|
||||
)
|
||||
for row in connection.execute(
|
||||
"SELECT device_id,link_id,measurement_unit,transmission_frequency "
|
||||
"FROM asset.scada_devices "
|
||||
"WHERE device_type IN ('pipe_flow','flow') AND device_id=ANY(%s)",
|
||||
([value.device_id for value in pipe.flow_bindings],),
|
||||
)
|
||||
}
|
||||
expected_pressure_bindings = {
|
||||
(
|
||||
value.device_id,
|
||||
value.node_id,
|
||||
expected.pressure_measurement_unit,
|
||||
str(timedelta(seconds=expected.report_timestep_seconds)),
|
||||
)
|
||||
for value in pipe.pressure_bindings
|
||||
}
|
||||
expected_flow_bindings = {
|
||||
(
|
||||
value.device_id,
|
||||
value.link_id,
|
||||
expected.flow_measurement_unit,
|
||||
str(timedelta(seconds=expected.report_timestep_seconds)),
|
||||
)
|
||||
for value in pipe.flow_bindings
|
||||
}
|
||||
if pressure_bindings != expected_pressure_bindings:
|
||||
raise ValueError(f"pipe {pipe.pipe_id} pressure bindings changed")
|
||||
if flow_bindings != expected_flow_bindings:
|
||||
raise ValueError(f"pipe {pipe.pipe_id} flow bindings changed")
|
||||
pressure_distances = {
|
||||
str(row["device_id"]): row["distance_m"]
|
||||
for row in connection.execute(
|
||||
"SELECT d.device_id,ST_Distance(d.geom,n.geom) AS distance_m "
|
||||
"FROM asset.scada_devices d "
|
||||
"JOIN gis.node_geometries n ON n.node_id=d.node_id "
|
||||
"WHERE d.device_id=ANY(%s)",
|
||||
([value.device_id for value in pipe.pressure_bindings],),
|
||||
)
|
||||
}
|
||||
flow_distances = {
|
||||
str(row["device_id"]): row["distance_m"]
|
||||
for row in connection.execute(
|
||||
"SELECT d.device_id,"
|
||||
"ST_Distance(d.geom,ST_Transform(p.geom,900914)) AS distance_m "
|
||||
"FROM asset.scada_devices d "
|
||||
"JOIN gis.pipes p ON p.id=d.link_id "
|
||||
"WHERE d.device_id=ANY(%s)",
|
||||
([value.device_id for value in pipe.flow_bindings],),
|
||||
)
|
||||
}
|
||||
if len(pressure_distances) != len(pipe.pressure_bindings) or any(
|
||||
distance is None
|
||||
or float(distance) > design.binding_validation.max_pressure_distance_m
|
||||
for distance in pressure_distances.values()
|
||||
):
|
||||
raise ValueError(f"pipe {pipe.pipe_id} pressure binding geometry changed")
|
||||
if len(flow_distances) != len(pipe.flow_bindings) or any(
|
||||
distance is None
|
||||
or float(distance) > design.binding_validation.max_flow_distance_m
|
||||
for distance in flow_distances.values()
|
||||
):
|
||||
raise ValueError(f"pipe {pipe.pipe_id} flow binding geometry changed")
|
||||
connection.rollback()
|
||||
|
||||
|
||||
def validate_tool_inputs(design) -> None:
|
||||
from app.mcp.extended_models import BurstSimulationInput
|
||||
|
||||
pipes = {pipe.pipe_id: pipe for pipe in design.pipes}
|
||||
for case in design.cases:
|
||||
pipe = pipes[case.pipe_id]
|
||||
for severity in design.calibration_order:
|
||||
BurstSimulationInput(
|
||||
context_handle="benchmark-design-context-handle-0001",
|
||||
timeout_seconds=240,
|
||||
start_time=case.occurred_at,
|
||||
duration=design.duration_seconds,
|
||||
output_node_ids=pipe.output_node_ids,
|
||||
output_link_ids=pipe.output_link_ids,
|
||||
burst_pipe_ids=[pipe.pipe_id],
|
||||
orifice_areas_cm2=[pipe.candidate_orifice_areas_cm2[severity]],
|
||||
)
|
||||
|
||||
|
||||
def main() -> int:
|
||||
parser = argparse.ArgumentParser(description=__doc__)
|
||||
parser.add_argument(
|
||||
"--design",
|
||||
type=Path,
|
||||
default=Path("resources/benchmarks/tjwater2-v2-burst-simulation.json"),
|
||||
)
|
||||
args = parser.parse_args()
|
||||
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
|
||||
|
||||
try:
|
||||
from scripts.benchmark.simulation_design import load_simulation_design
|
||||
|
||||
design = load_simulation_design(args.design)
|
||||
validate_tool_inputs(design)
|
||||
biz_dsn, _ = asyncio.run(resolve_project(design))
|
||||
physical_name = conninfo_to_dict(biz_dsn).get("dbname")
|
||||
if not physical_name:
|
||||
raise ValueError("business route has no database name")
|
||||
template_dsn = make_conninfo(biz_dsn, dbname=physical_name + "_template")
|
||||
with readonly(biz_dsn) as connection:
|
||||
validate_database(connection, design, template=False)
|
||||
with readonly(template_dsn) as connection:
|
||||
validate_database(connection, design, template=True)
|
||||
except Exception as exc:
|
||||
frames = traceback.extract_tb(exc.__traceback__)
|
||||
location = frames[-1] if frames else None
|
||||
suffix = f" at {location.filename}:{location.lineno}" if location else ""
|
||||
print(f"burst simulation design check failed ({type(exc).__name__}){suffix}", file=sys.stderr)
|
||||
return 1
|
||||
|
||||
maximum_runs = len(design.cases) * len(design.severity_levels)
|
||||
print(
|
||||
f"burst simulation design ready: {len(design.pipes)} pipes, "
|
||||
f"{len(design.cases)} cases, at most {maximum_runs} calibration runs; "
|
||||
"tool inputs, SCADA bindings, BizDB and model template were read-only verified"
|
||||
)
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main())
|
||||
Reference in New Issue
Block a user