From 5c6ed694dadb06b35b05ce57eb7d9bdcad03b4ab Mon Sep 17 00:00:00 2001 From: Dan Martin <93075337+Mingqwqqaq@users.noreply.github.com> Date: Wed, 19 Aug 2026 12:44:37 +0800 Subject: [PATCH 1/2] feat(WastewaterTreatment): add BSM1 aeration control benchmark --- TASK_DETAILS.md | 5 + TASK_DETAILS_zh-CN.md | 5 + .../BSM1AerationControl/README.md | 53 +++ .../BSM1AerationControl/README_zh-CN.md | 44 +++ .../BSM1AerationControl/Task.md | 92 +++++ .../BSM1AerationControl/Task_zh-CN.md | 58 +++ .../baseline/result_log.txt | 24 ++ .../BSM1AerationControl/baseline/solution.py | 16 + .../frontier_eval/agent_files.txt | 7 + .../frontier_eval/artifact_files.txt | 2 + .../frontier_eval/candidate_destination.txt | 1 + .../frontier_eval/constraints.txt | 9 + .../frontier_eval/copy_files.txt | 12 + .../frontier_eval/eval_command.txt | 1 + .../frontier_eval/initial_program.txt | 1 + .../frontier_eval/readonly_files.txt | 15 + .../references/BSD-3-Clause-bsm2-python.txt | 28 ++ .../references/config.json | 63 ++++ .../references/design_notes.md | 41 +++ .../BSM1AerationControl/scripts/init.py | 20 + .../verification/bsm1_model.py | 341 ++++++++++++++++++ .../verification/evaluator.py | 309 ++++++++++++++++ .../verification/policy_runtime.py | 251 +++++++++++++ .../verification/policy_worker.py | 98 +++++ .../verification/requirements.txt | 2 + .../verification/test_bsm1_model.py | 50 +++ .../verification/test_evaluator.py | 81 +++++ .../verification/test_policy_runtime.py | 46 +++ benchmarks/WastewaterTreatment/README.md | 11 + .../WastewaterTreatment/README_zh-CN.md | 9 + .../conf/task/bsm1_aeration_control.yaml | 6 + 31 files changed, 1701 insertions(+) create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/README.md create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/README_zh-CN.md create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/Task.md create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/Task_zh-CN.md create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/solution.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/agent_files.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/artifact_files.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/candidate_destination.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/constraints.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/copy_files.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/eval_command.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/initial_program.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/readonly_files.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/references/BSD-3-Clause-bsm2-python.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/references/config.json create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/scripts/init.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/bsm1_model.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/evaluator.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_runtime.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/requirements.txt create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_evaluator.py create mode 100644 benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py create mode 100644 benchmarks/WastewaterTreatment/README.md create mode 100644 benchmarks/WastewaterTreatment/README_zh-CN.md create mode 100644 frontier_eval/conf/task/bsm1_aeration_control.yaml diff --git a/TASK_DETAILS.md b/TASK_DETAILS.md index 7475c4cc..2ff811fb 100644 --- a/TASK_DETAILS.md +++ b/TASK_DETAILS.md @@ -351,5 +351,10 @@ We welcome new engineering problem ideas — even without complete verification DiffSimThermalControl Process optimization in additive manufacturing via differentiable simulation + + WastewaterTreatment + BSM1AerationControl + Feedback control of activated-sludge aeration and internal recycle across dry, rain, and storm operation + diff --git a/TASK_DETAILS_zh-CN.md b/TASK_DETAILS_zh-CN.md index e2a070a2..e3a4fa6c 100644 --- a/TASK_DETAILS_zh-CN.md +++ b/TASK_DETAILS_zh-CN.md @@ -351,5 +351,10 @@ Frontier-Eng 目前已覆盖以下领域的任务。每个任务均配有可运 DiffSimThermalControl 基于可微仿真的增材制造工艺优化 + + WastewaterTreatment + BSM1AerationControl + 在旱天、降雨和暴雨工况下反馈控制活性污泥曝气与内回流 + diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/README.md b/benchmarks/WastewaterTreatment/BSM1AerationControl/README.md new file mode 100644 index 00000000..aeca69ee --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/README.md @@ -0,0 +1,53 @@ +# BSM1 Aeration Control + +Design a deterministic feedback controller for an activated-sludge plant derived from the +IWA Benchmark Simulation Model No. 1 (BSM1). The controller sets oxygen-transfer coefficients +in the three aerobic reactors and the internal recycle flow. It must balance effluent quality, +energy, compliance, and actuator smoothness across dry, rain, and storm scenarios. + +Edit only the EVOLVE-BLOCK in `scripts/init.py`, preserving: + +```python +def reset_controller(scenario: dict) -> None: ... +def control(observation: dict) -> dict: ... +``` + +## Setup + +The task is offline and CPU-only: + +```bash +python -m pip install -r verification/requirements.txt +``` + +No raw IWA influent files are redistributed. The evaluator generates deterministic trajectories +from published BSM1 averages and weather-event descriptions. A complete baseline run takes about +45 seconds on a laptop. + +## Direct evaluation + +```bash +python verification/evaluator.py scripts/init.py --metrics-out metrics.json --artifacts-out artifacts.json +``` + +## Regression tests + +```bash +python -m unittest discover -s verification -p "test_*.py" -v +``` + +## Unified evaluation + +From the repository root: + +```bash +python -m frontier_eval task=unified task.benchmark=WastewaterTreatment/BSM1AerationControl algorithm=openevolve algorithm.iterations=0 +``` + +The ranking metric is `combined_score` (higher is better). `metrics.json` contains the aggregate +score and validity fields; `artifacts.json` contains scenario-level engineering metrics and daily +effluent samples. Candidate code runs in a separate, bounded JSON-lines worker process. This is +process isolation, not an operating-system security sandbox. + +See `Task.md` for the exact interface and scoring model, and +`references/design_notes.md` for provenance, validation, and modelling limitations. diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/README_zh-CN.md b/benchmarks/WastewaterTreatment/BSM1AerationControl/README_zh-CN.md new file mode 100644 index 00000000..baa70329 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/README_zh-CN.md @@ -0,0 +1,44 @@ +# BSM1 曝气控制 + +为基于 IWA 一号基准仿真模型（BSM1）的活性污泥处理厂设计确定性反馈控制器。控制器设置三个好氧反应池的氧传质系数与内回流量，并在旱天、降雨和暴雨场景下权衡出水质量、能耗、达标情况和执行器平滑性。 + +只修改 `scripts/init.py` 中的 EVOLVE-BLOCK，并保持以下接口： + +```python +def reset_controller(scenario: dict) -> None: ... +def control(observation: dict) -> dict: ... +``` + +## 环境安装 + +本任务离线运行且只需要 CPU： + +```bash +python -m pip install -r verification/requirements.txt +``` + +仓库不再分发 IWA 原始进水文件。评测器依据已发表的 BSM1 平均值与天气事件说明生成确定性轨迹。在普通笔记本上完整基准评测约需 45 秒。 + +## 直接评测 + +```bash +python verification/evaluator.py scripts/init.py --metrics-out metrics.json --artifacts-out artifacts.json +``` + +## 回归测试 + +```bash +python -m unittest discover -s verification -p "test_*.py" -v +``` + +## 统一评测 + +在仓库根目录运行： + +```bash +python -m frontier_eval task=unified task.benchmark=WastewaterTreatment/BSM1AerationControl algorithm=openevolve algorithm.iterations=0 +``` + +排名指标为 `combined_score`，越高越好。`metrics.json` 给出聚合分数与有效性；`artifacts.json` 给出逐场景工程指标和每日出水样本。候选代码运行在带时限的独立 JSON-lines 工作进程中；这是进程隔离，不是操作系统级安全沙箱。 + +完整接口和评分模型见 `Task_zh-CN.md`，来源、验证和模型限制见 `references/design_notes.md`。 diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/Task.md b/benchmarks/WastewaterTreatment/BSM1AerationControl/Task.md new file mode 100644 index 00000000..0e2e2ebf --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/Task.md @@ -0,0 +1,92 @@ +# Task: BSM1 Aeration Control + +## Engineering setting + +Activated-sludge plants must remove carbon and nitrogen while limiting blower and pumping energy. +Influent flow and composition vary diurnally and become more difficult during rain and storm +events. This task exposes feedback control rather than a one-shot parameter fit: every 15 minutes, +the policy observes noisy process measurements and selects aeration and internal recycle settings. + +## Process model and scenarios + +The evaluator implements the standard BSM1 layout: five completely mixed ASM1 reactors (two +anoxic, three aerobic) followed by a ten-layer Takacs secondary settler. Fixed return-activated +sludge and waste-sludge flows are 18,446 and 385 m3/day. The plant is simulated for 14 days and +the final seven days are scored. + +Three deterministic scenarios share the same diurnal base load: + +- `dry`: diurnal and weekend load variation; +- `rain`: a sustained dilution-water event from day 8.35 to day 10.44; +- `storm`: two shorter hydraulic pulses centred near days 8.87 and 11.18. + +The rain and storm series are deterministic IWA-derived engineering trajectories, not byte-for-byte +copies of the official BSM1 dynamic influent files. Sensor noise uses fixed, scenario-specific seeds. + +## Observation + +`control(observation)` receives only current or past information: + +- `scenario_id`, `weather`, `time_day`, and `step_minutes`; +- influent flow and ammonium; +- dissolved oxygen in reactors 3, 4, and 5; +- nitrate in reactor 2; +- effluent ammonium and total nitrogen; +- `previous_action`. + +The dissolved-oxygen, nitrate, and ammonium measurements contain deterministic sensor noise. The +policy does not receive future influent, process state arrays, random seeds, or evaluator internals. + +## Action interface and hard constraints + +Return exactly these four finite numeric fields: + +```python +{ + "kla3_per_day": float, # [0, 360] + "kla4_per_day": float, # [0, 360] + "kla5_per_day": float, # [0, 360] + "internal_recycle_m3_per_day": float, # [0, 92230] +} +``` + +At one 15-minute step, each KLa may change by at most 120/day and internal recycle by at most +30,000 m3/day. Missing/extra fields, booleans, non-finite values, range errors, slew errors, +exceptions, protocol errors, or timeouts make the complete candidate invalid with score zero. +`reset_controller(scenario)` must reset all candidate-owned state before each scenario. + +## Metrics and score + +For every scored step, the evaluator computes the official BSM1 effluent quality index (EQI), +aeration energy, pumping energy, mixing energy, actuator switching, and normalized exceedance of +five standard limits: NH4-N 4, total nitrogen 18, COD 100, TSS 30, and BOD5 10 g/m3. + +Each objective is mapped to a dimensionless utility: + +```text +Uq = exp(-EQI / 6000) +Ua = exp(-aeration_energy / 5000) +Up = exp(-pumping_energy / 1200) +Um = exp(-mixing_energy / 600) +Us = exp(-4 * switching_index) + +base = 100 * (0.45 Uq + 0.20 Ua + 0.10 Up + 0.10 Um + 0.15 Us) +scenario_score = base * exp(-6 * violation_index) +``` + +The aggregate rewards typical and worst-case behavior: + +```text +combined_score = 0.75 * mean(scenario_scores) + 0.25 * min(scenario_scores) +``` + +The score is absolute and contains no frozen candidate baseline. A fixed published BSM1 operating +point is supplied only as editable starter code and as a reproducible comparison in the result log. + +## Candidate process boundary + +The evaluator copies the candidate to a temporary directory and imports it in a persistent worker. +Communication is JSON-only; import, each call, cumulative response time, and response size are +bounded. Candidate stdout is discarded and credentials are removed from its environment. This +prevents ordinary Python-level mutation of the parent evaluator; it is not a sandbox against +hostile native code or unrestricted filesystem/network access. diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/Task_zh-CN.md b/benchmarks/WastewaterTreatment/BSM1AerationControl/Task_zh-CN.md new file mode 100644 index 00000000..a0720228 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/Task_zh-CN.md @@ -0,0 +1,58 @@ +# 任务：BSM1 曝气控制 + +## 工程背景 + +活性污泥处理厂既要去除有机物和氮，又要控制鼓风与泵送能耗。进水流量和组分具有日周期变化，降雨与暴雨会进一步增加水力冲击。本任务是反馈控制而非一次性参数拟合：策略每 15 分钟读取带噪过程测量，并设置曝气与内回流。 + +## 工艺模型与场景 + +评测器实现标准 BSM1 布局：五个完全混合 ASM1 反应池（前两个缺氧、后三个好氧）以及十层 Takacs 二沉池。固定污泥回流和排泥流量分别为 18,446 与 385 m3/day。每个场景仿真 14 天，最后 7 天计分。 + +三个确定性场景共享同一日周期基础负荷： + +- `dry`：日周期和周末负荷变化； +- `rain`：第 8.35 至 10.44 天的持续稀释水事件； +- `storm`：中心约在第 8.87 与 11.18 天的两次短时水力脉冲。 + +降雨与暴雨序列是依据 IWA 资料构造的确定性工程轨迹，并非官方 BSM1 动态进水文件的逐字节副本。传感器噪声使用固定且按场景区分的种子。 + +## 观测 + +`control(observation)` 只接收当前或历史信息：场景、天气、时间、步长；进水流量与氨氮；3—5 号池溶解氧；2 号池硝酸盐；出水氨氮与总氮；以及 `previous_action`。溶解氧、硝酸盐和氨氮测量含确定性噪声。策略不能获得未来进水、完整过程状态、随机种子或评测器内部信息。 + +## 动作接口与硬约束 + +必须恰好返回四个有限数值字段： + +```python +{ + "kla3_per_day": float, # [0, 360] + "kla4_per_day": float, # [0, 360] + "kla5_per_day": float, # [0, 360] + "internal_recycle_m3_per_day": float, # [0, 92230] +} +``` + +相邻 15 分钟内，每个 KLa 的最大变化为 120/day，内回流最大变化为 30,000 m3/day。字段缺失或多余、布尔值、非有限值、越界、变化率超限、异常、协议错误或超时都会使整个候选无效并得零分。`reset_controller(scenario)` 必须在每个场景前重置策略自身状态。 + +## 指标与评分 + +计分期的每一步都会计算 BSM1 出水质量指数（EQI）、曝气/泵送/搅拌能耗、执行器切换，以及五项标准限值的归一化超标量：NH4-N 4、总氮 18、COD 100、TSS 30、BOD5 10 g/m3。 + +```text +Uq = exp(-EQI / 6000) +Ua = exp(-曝气能耗 / 5000) +Up = exp(-泵送能耗 / 1200) +Um = exp(-搅拌能耗 / 600) +Us = exp(-4 * switching_index) + +base = 100 * (0.45 Uq + 0.20 Ua + 0.10 Up + 0.10 Um + 0.15 Us) +场景分 = base * exp(-6 * violation_index) +combined_score = 0.75 * 场景均分 + 0.25 * 最差场景分 +``` + +该分数是绝对工程分，不依赖冻结的候选基线。公开的 BSM1 固定工况仅作为可编辑起始程序和结果日志中的可复现实验对照。 + +## 候选进程边界 + +评测器把候选复制到临时目录并在持久工作进程中导入。通信仅使用 JSON；导入、单次调用、累计响应时间和响应大小均受限。候选标准输出会被丢弃，环境中的凭据会被移除。该边界可阻止普通 Python 代码修改父评测器，但不能替代针对恶意原生代码、文件系统或网络访问的操作系统沙箱。 diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt new file mode 100644 index 00000000..954cc0da --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt @@ -0,0 +1,24 @@ +Direct command: +python verification/evaluator.py scripts/init.py --metrics-out metrics.json --artifacts-out artifacts.json + +Environment: +Windows 11, Python 3.11, CPU-only + +Result: +=== BSM1 Aeration Control === +scenario=dry score=53.457 eqi=5488.86 ae=3341.39 violations=0.011330 +scenario=rain score=44.609 eqi=7535.45 ae=3341.39 violations=0.025575 +scenario=storm score=48.513 eqi=6780.93 ae=3341.39 violations=0.017011 +--- +completed_scenarios: 3/3 +diagnostic_score: 47.7973 +combined_score: 47.7973 + +Regression: +python -m unittest discover -s verification -p "test_*.py" -v +Ran 9 tests in 3.574s +OK + +Unified evaluation: +python -m frontier_eval task=unified task.benchmark=WastewaterTreatment/BSM1AerationControl algorithm=openevolve algorithm.iterations=0 +combined_score=47.7973, valid=1.0000, completed_scenarios=3.0000, benchmark_returncode=0.0000 diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/solution.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/solution.py new file mode 100644 index 00000000..5911e71a --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/solution.py @@ -0,0 +1,16 @@ +"""Reference fixed-operating-point policy for BSM1AerationControl.""" + +from __future__ import annotations + + +def reset_controller(scenario: dict) -> None: + pass + + +def control(observation: dict) -> dict: + return { + "kla3_per_day": 240.0, + "kla4_per_day": 240.0, + "kla5_per_day": 84.0, + "internal_recycle_m3_per_day": 55338.0, + } diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/agent_files.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/agent_files.txt new file mode 100644 index 00000000..e8a78865 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/agent_files.txt @@ -0,0 +1,7 @@ +README.md +README_zh-CN.md +Task.md +Task_zh-CN.md +scripts/init.py +references/design_notes.md +frontier_eval/constraints.txt diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/artifact_files.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/artifact_files.txt new file mode 100644 index 00000000..76dc893a --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/artifact_files.txt @@ -0,0 +1,2 @@ +metrics.json +artifacts.json diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/candidate_destination.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/candidate_destination.txt new file mode 100644 index 00000000..b9411b3d --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/candidate_destination.txt @@ -0,0 +1 @@ +scripts/init.py diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/constraints.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/constraints.txt new file mode 100644 index 00000000..1b10a473 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/constraints.txt @@ -0,0 +1,9 @@ +BSM1 aeration-control constraints: +1) Edit only `scripts/init.py` inside the EVOLVE-BLOCK markers. +2) Keep `reset_controller(scenario) -> None` and `control(observation) -> dict` working. +3) Return exactly the four finite numeric fields documented in Task.md. +4) Keep each KLa in [0, 360]/day and internal recycle in [0, 92230] m3/day. +5) Respect per-step slew limits: 120/day for each KLa and 30000 m3/day for recycle. +6) Reset all candidate-owned state before each dry, rain, or storm scenario. +7) Keep the controller deterministic, self-contained, and within the 0.10-second call limit. +8) Do not read, write, import, or modify evaluator, model, result, or environment-secret files. diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/copy_files.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/copy_files.txt new file mode 100644 index 00000000..68c8d26c --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/copy_files.txt @@ -0,0 +1,12 @@ +README.md +README_zh-CN.md +Task.md +Task_zh-CN.md +scripts/init.py +verification/evaluator.py +verification/bsm1_model.py +verification/policy_runtime.py +verification/policy_worker.py +verification/requirements.txt +references/config.json +references/BSD-3-Clause-bsm2-python.txt diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/eval_command.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/eval_command.txt new file mode 100644 index 00000000..527d23b5 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/eval_command.txt @@ -0,0 +1 @@ +{python} verification/evaluator.py {candidate} --metrics-out metrics.json --artifacts-out artifacts.json diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/initial_program.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/initial_program.txt new file mode 100644 index 00000000..b9411b3d --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/initial_program.txt @@ -0,0 +1 @@ +scripts/init.py diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/readonly_files.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/readonly_files.txt new file mode 100644 index 00000000..ce96ab94 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/frontier_eval/readonly_files.txt @@ -0,0 +1,15 @@ +README.md +README_zh-CN.md +Task.md +Task_zh-CN.md +verification/evaluator.py +verification/bsm1_model.py +verification/policy_runtime.py +verification/policy_worker.py +verification/test_bsm1_model.py +verification/test_evaluator.py +verification/test_policy_runtime.py +verification/requirements.txt +references/config.json +references/design_notes.md +references/BSD-3-Clause-bsm2-python.txt diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/references/BSD-3-Clause-bsm2-python.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/BSD-3-Clause-bsm2-python.txt new file mode 100644 index 00000000..eb25b478 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/BSD-3-Clause-bsm2-python.txt @@ -0,0 +1,28 @@ +BSD 3-Clause License + +Copyright (c) 2025, FAU-EVT + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/references/config.json b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/config.json new file mode 100644 index 00000000..fe7018fa --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/config.json @@ -0,0 +1,63 @@ +{ + "benchmark": "BSM1AerationControl", + "source": { + "iwa_page": "https://iwa-mia.org/benchmarking/", + "iwa_report": "https://iwa-mia.org/wp-content/uploads/2019/04/BSM_TG_Tech_Report_no_1_BSM1_General_Description.pdf", + "validation_implementation": "https://github.com/fau-evt/bsm2-python", + "validation_license": "BSD-3-Clause" + }, + "simulation": { + "days": 14.0, + "evaluation_start_day": 7.0, + "step_minutes": 15.0, + "sensor_noise_seed": 1729, + "sensor_noise_std": { + "dissolved_oxygen": 0.025, + "nitrate": 0.025, + "ammonium": 0.04 + } + }, + "scenarios": [ + {"scenario_id": "dry", "weather": "dry"}, + {"scenario_id": "rain", "weather": "rain"}, + {"scenario_id": "storm", "weather": "storm"} + ], + "action": { + "kla_min_per_day": 0.0, + "kla_max_per_day": 360.0, + "qintr_min_m3_per_day": 0.0, + "qintr_max_m3_per_day": 92230.0, + "kla_max_change_per_step": 120.0, + "qintr_max_change_per_step": 30000.0 + }, + "reference_action": { + "kla3_per_day": 240.0, + "kla4_per_day": 240.0, + "kla5_per_day": 84.0, + "internal_recycle_m3_per_day": 55338.0 + }, + "limits": { + "ammonium_gN_per_m3": 4.0, + "total_nitrogen_gN_per_m3": 18.0, + "cod_gCOD_per_m3": 100.0, + "tss_gSS_per_m3": 30.0, + "bod5_gBOD_per_m3": 10.0 + }, + "scoring": { + "eqi_scale": 6000.0, + "aeration_energy_scale": 5000.0, + "pumping_energy_scale": 1200.0, + "mixing_energy_scale": 600.0, + "switching_exponent": 4.0, + "violation_exponent": 6.0, + "objective_weights": { + "quality": 0.45, + "aeration": 0.20, + "pumping": 0.10, + "mixing": 0.10, + "smoothness": 0.15 + }, + "mean_weight": 0.75, + "worst_case_weight": 0.25 + } +} diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md new file mode 100644 index 00000000..10855221 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md @@ -0,0 +1,41 @@ +# Model provenance and design notes + +## Sources + +The process layout, ASM1 kinetics, Takacs settler, standard operating point, effluent limits, +effluent quality index, and energy equations follow the IWA Benchmark Simulation Model No. 1: + +- [IWA benchmark overview](https://iwa-mia.org/benchmarking/) +- [BSM1 general description and simulation protocol](https://iwa-mia.org/wp-content/uploads/2019/04/BSM_TG_Tech_Report_no_1_BSM1_General_Description.pdf) + +The numerical implementation was cross-checked against the open-source +[`fau-evt/bsm2-python`](https://github.com/fau-evt/bsm2-python) implementation. Initial +steady-state arrays and equation structure were adapted from that BSD-3-Clause project; its +license is retained in `BSD-3-Clause-bsm2-python.txt`. + +## Data boundary + +The official IWA site publishes dynamic influent files but does not state a redistribution license +for those data. This repository therefore does not copy them. `influent_at()` builds reproducible +engineering trajectories from the published average BSM1 influent, a deterministic diurnal/weekend +profile, and the dry/rain/storm timing described by IWA. Rain and storm add water while conserving +the instantaneous mass flow of all 13 ASM1 components. Results should be described as BSM1-derived, +not as an exact replay of the official influent files. + +## Validation evidence + +- The fixed operating point yields 3341.3867 kWh/day aeration energy, matching the BSM1 reference + equation and published operating point. +- Unit tests check weather-event component mass conservation, positive settler hydraulics after a + recycle-flow change, action validation, deterministic results, worker isolation, and timeouts. +- A 96-step dry-weather regression demonstrates optimization headroom: lowering internal recycle + reduces pumping energy and improves score relative to the starter policy. +- The full direct evaluation covers all three 14-day scenarios; its recorded output is in + `baseline/result_log.txt`. + +## Deliberate limitations + +This compact CPU benchmark assumes ideal completely mixed reactors, a fixed temperature of 15 C, +fixed return/waste sludge flows, no explicit blower or pump dynamics, and deterministic hydraulic +events. Candidate actions cannot change exogenous influent or sensor-noise sequences. The model is +appropriate for controller-search benchmarking, not plant commissioning or regulatory prediction. diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/scripts/init.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/scripts/init.py new file mode 100644 index 00000000..e71ecf02 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/scripts/init.py @@ -0,0 +1,20 @@ +"""Editable baseline controller for the BSM1 aeration benchmark.""" + +from __future__ import annotations + + +def reset_controller(scenario: dict) -> None: + """Reset any candidate-owned state before a weather scenario.""" + + +def control(observation: dict) -> dict: + """Return aeration coefficients and the internal recycle flow.""" + + # EVOLVE-BLOCK-START + return { + "kla3_per_day": 240.0, + "kla4_per_day": 240.0, + "kla5_per_day": 84.0, + "internal_recycle_m3_per_day": 55338.0, + } + # EVOLVE-BLOCK-END diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/bsm1_model.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/bsm1_model.py new file mode 100644 index 00000000..bb4ea899 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/bsm1_model.py @@ -0,0 +1,341 @@ +"""Independent BSM1 plant model used by the aeration-control benchmark. + +The equations follow the IWA Benchmark Simulation Model No. 1: five ASM1 +reactors and a ten-layer Takacs secondary settler. The implementation was +cross-checked against ``fau-evt/bsm2-python`` and is adapted from its +BSD-3-Clause implementation; the retained license is in ``references/``. + +Units are days, m3, g/m3, and kWh/day unless stated otherwise. +""" + +from __future__ import annotations + +import json +import math +from functools import lru_cache +from pathlib import Path +from typing import Any, Mapping + +import numpy as np +from scipy.integrate import odeint + + +ROOT = Path(__file__).resolve().parents[1] +CONFIG_PATH = ROOT / "references" / "config.json" + +SI, SS, XI, XS, XBH, XBA, XP, SO, SNO, SNH, SND, XND, SALK = range(13) +TSS, Q, TEMP = 13, 14, 15 +N_COMPONENTS = 16 + +ASM1_PARAMETERS = np.array( + [ + 4.0, 10.0, 0.2, 0.5, 0.3, 0.5, 1.0, 0.4, 0.05, 0.8, + 0.05, 3.0, 0.1, 0.8, 0.67, 0.24, 0.08, 0.08, 0.06, + 0.75, 0.75, 0.75, 0.75, 0.75, + ], + dtype=float, +) +VOLUMES = np.array([1000.0, 1000.0, 1333.0, 1333.0, 1333.0]) +OXYGEN_SATURATION = 8.0 +RETURN_FLOW = 18446.0 +WASTE_FLOW = 385.0 +SETTLER_AREA = 1500.0 +SETTLER_HEIGHT = 4.0 +SETTLER_LAYERS = 10 +SETTLER_FEED_LAYER = 5 +SETTLER_PARAMETERS = np.array([250.0, 474.0, 0.000576, 0.00286, 0.00228, 3000.0, 3000.0]) + +_REACTOR_INITIAL = np.array( + [ + [30.0, 2.696169337, 1149.16340146, 79.78436831, 2553.45921988, 151.75632721, 449.22675116, 0.018822826, 8.310656768, 7.108011983, 1.228022799, 5.127002559, 4.659811087], + [30.0, 1.402154895, 1149.16340144, 73.67182578, 2555.28394130, 151.67460865, 449.89826272, 0.000279954, 6.570406004, 7.536848342, 0.896501365, 4.847841192, 4.814745881], + [30.0, 1.113091672, 1149.16340141, 62.63753037, 2558.68270986, 152.32053441, 450.79460391, 2.0, 9.544077377, 4.687308245, 0.822255969, 4.244269898, 4.398802205], + [30.0, 0.967558516, 1149.16340138, 54.01878920, 2560.36433241, 152.85355533, 451.69155922, 2.0, 12.068200287, 2.322898686, 0.751729972, 3.767047773, 4.049621314], + [30.0, 0.857137129, 1149.16340135, 47.36389426, 2560.70622587, 153.17106513, 452.58865148, 2.0, 13.711131912, 0.866431079, 0.690377137, 3.394128545, 3.828235655], + ], + dtype=float, +) +_SETTLER_TSS_INITIAL = np.array( + [12.50109702, 18.11774069, 29.54712298, 68.99884217, 356.2589969, + 356.25899691, 356.2589969, 356.25899691, 356.2589969, 6398.69776357], + dtype=float, +) + + +@lru_cache(maxsize=1) +def load_config() -> dict[str, Any]: + return json.loads(CONFIG_PATH.read_text(encoding="utf-8")) + + +def _tss(components: np.ndarray) -> float: + return float(0.75 * np.sum(components[[XI, XS, XBH, XBA, XP]])) + + +def _stream(components: np.ndarray, flow: float, temperature: float = 15.0) -> np.ndarray: + result = np.zeros(N_COMPONENTS, dtype=float) + result[:13] = components + result[TSS] = _tss(components) + result[Q] = float(flow) + result[TEMP] = float(temperature) + return result + + +def combine_streams(*streams: np.ndarray) -> np.ndarray: + total_flow = float(sum(max(0.0, stream[Q]) for stream in streams)) + if total_flow <= 0.0: + raise ValueError("cannot combine zero-flow streams") + result = np.zeros(N_COMPONENTS, dtype=float) + for stream in streams: + result[:14] += stream[:14] * stream[Q] + result[TEMP] += stream[TEMP] * stream[Q] + result[:14] /= total_flow + result[TEMP] /= total_flow + result[Q] = total_flow + return result + + +def _asm1_rhs(y: np.ndarray, influent: np.ndarray, kla: float, volume: float) -> np.ndarray: + y = np.maximum(np.asarray(y, dtype=float), 1e-12) + p = ASM1_PARAMETERS + mu_h, k_s, k_oh, k_no, b_h, mu_a, k_nh, k_oa, b_a = p[:9] + ny_g, k_a, k_h, k_x, ny_h, y_h, y_a, f_p, i_xb, i_xp = p[9:19] + temp = float(influent[TEMP]) + mu_h *= math.exp((math.log(mu_h / 3.0) / 5.0) * (temp - 15.0)) + b_h *= math.exp((math.log(b_h / 0.2) / 5.0) * (temp - 15.0)) + mu_a *= math.exp((math.log(mu_a / 0.3) / 5.0) * (temp - 15.0)) + b_a *= math.exp((math.log(b_a / 0.03) / 5.0) * (temp - 15.0)) + k_h *= math.exp((math.log(k_h / 2.5) / 5.0) * (temp - 15.0)) + k_a *= math.exp((math.log(k_a / 0.04) / 5.0) * (temp - 15.0)) + kla_temp = kla * 1.024 ** (temp - 15.0) + + proc1 = mu_h * y[SS] / (k_s + y[SS]) * y[SO] / (k_oh + y[SO]) * y[XBH] + proc2 = mu_h * y[SS] / (k_s + y[SS]) * k_oh / (k_oh + y[SO]) * y[SNO] / (k_no + y[SNO]) * ny_g * y[XBH] + proc3 = mu_a * y[SNH] / (k_nh + y[SNH]) * y[SO] / (k_oa + y[SO]) * y[XBA] + proc4 = b_h * y[XBH] + proc5 = b_a * y[XBA] + proc6 = k_a * y[SND] * y[XBH] + ratio_xs = y[XS] / max(y[XBH], 1e-12) + proc7 = k_h * ratio_xs / (k_x + ratio_xs) * ( + y[SO] / (k_oh + y[SO]) + ny_h * k_oh / (k_oh + y[SO]) * y[SNO] / (k_no + y[SNO]) + ) * y[XBH] + proc8 = proc7 * y[XND] / max(y[XS], 1e-12) + + reaction = np.zeros(13, dtype=float) + reaction[SS] = (-proc1 - proc2) / y_h + proc7 + reaction[XS] = (1.0 - f_p) * (proc4 + proc5) - proc7 + reaction[XBH] = proc1 + proc2 - proc4 + reaction[XBA] = proc3 - proc5 + reaction[XP] = f_p * (proc4 + proc5) + reaction[SO] = -(1.0 - y_h) / y_h * proc1 - (4.57 - y_a) / y_a * proc3 + reaction[SNO] = -(1.0 - y_h) / (2.86 * y_h) * proc2 + proc3 / y_a + reaction[SNH] = -i_xb * (proc1 + proc2) - (i_xb + 1.0 / y_a) * proc3 + proc6 + reaction[SND] = -proc6 + proc8 + reaction[XND] = (i_xb - f_p * i_xp) * (proc4 + proc5) - proc8 + reaction[SALK] = ( + -i_xb / 14.0 * proc1 + + ((1.0 - y_h) / (14.0 * 2.86 * y_h) - i_xb / 14.0) * proc2 + - (i_xb / 14.0 + 1.0 / (7.0 * y_a)) * proc3 + + proc6 / 14.0 + ) + derivative = influent[Q] / volume * (influent[:13] - y) + reaction + derivative[SO] += kla_temp * (OXYGEN_SATURATION - y[SO]) + return derivative + + +def _integrate_reactor(state: np.ndarray, influent: np.ndarray, kla: float, volume: float, dt: float) -> np.ndarray: + result = odeint( + lambda y, _: _asm1_rhs(y, influent, kla, volume), + state, + np.array([0.0, dt]), + rtol=1e-4, + atol=1e-6, + mxstep=500, + )[-1] + if not np.all(np.isfinite(result)): + raise FloatingPointError("ASM1 integration produced non-finite state") + return np.maximum(result, 0.0) + + +def _settler_rhs(state: np.ndarray, feed: np.ndarray) -> np.ndarray: + layers = SETTLER_LAYERS + feed_layer = SETTLER_FEED_LAYER + height = SETTLER_HEIGHT / layers + state = np.maximum(np.asarray(state, dtype=float), 1e-8) + q_under = RETURN_FLOW + WASTE_FLOW + q_effluent = feed[Q] - q_under + if q_effluent <= 0.0: + raise ValueError("settler effluent flow is non-positive") + v_in = feed[Q] / SETTLER_AREA + v_up = q_effluent / SETTLER_AREA + v_down = q_under / SETTLER_AREA + + tss = state[7 * layers : 8 * layers] + v0_max, v0, r_h, r_p, f_ns, x_threshold, _ = SETTLER_PARAMETERS + velocity = v0 * ( + np.exp(-r_h * (tss - f_ns * feed[TSS])) + - np.exp(-r_p * (tss - f_ns * feed[TSS])) + ) + velocity = np.clip(velocity, 0.0, v0_max) + raw_flux = velocity * tss + flux = np.zeros(layers + 1, dtype=float) + for index in range(layers - 1): + if index < feed_layer - 1 and tss[index + 1] <= x_threshold: + flux[index + 1] = raw_flux[index] + else: + flux[index + 1] = min(raw_flux[index], raw_flux[index + 1]) + + derivative = np.zeros_like(state) + soluble_indices = [SI, SS, SO, SNO, SNH, SND, SALK] + for block, component in enumerate(soluble_indices): + values = state[block * layers : (block + 1) * layers] + out = derivative[block * layers : (block + 1) * layers] + for index in range(feed_layer - 1): + out[index] = v_up * (values[index + 1] - values[index]) / height + f = feed_layer - 1 + out[f] = (v_in * feed[component] - (v_up + v_down) * values[f]) / height + for index in range(feed_layer, layers): + out[index] = v_down * (values[index - 1] - values[index]) / height + + out_tss = derivative[7 * layers : 8 * layers] + for index in range(feed_layer - 1): + out_tss[index] = (v_up * (tss[index + 1] - tss[index]) - flux[index + 1] + flux[index]) / height + f = feed_layer - 1 + out_tss[f] = (v_in * feed[TSS] - (v_up + v_down) * tss[f] - flux[f + 1] + flux[f]) / height + for index in range(feed_layer, layers): + out_tss[index] = (v_down * (tss[index - 1] - tss[index]) - flux[index + 1] + flux[index]) / height + return derivative + + +def _settler_outputs(state: np.ndarray, feed: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + layers = SETTLER_LAYERS + state = np.maximum(state, 0.0) + soluble_indices = [SI, SS, SO, SNO, SNH, SND, SALK] + effluent = np.zeros(13, dtype=float) + return_sludge = np.zeros(13, dtype=float) + for block, component in enumerate(soluble_indices): + values = state[block * layers : (block + 1) * layers] + effluent[component] = values[0] + return_sludge[component] = values[-1] + top_tss = state[7 * layers] + bottom_tss = state[8 * layers - 1] + particulate = [XI, XS, XBH, XBA, XP, XND] + if feed[TSS] > 1e-12: + for component in particulate: + effluent[component] = top_tss / feed[TSS] * feed[component] + return_sludge[component] = bottom_tss / feed[TSS] * feed[component] + return _stream(effluent, feed[Q] - RETURN_FLOW - WASTE_FLOW), _stream(return_sludge, RETURN_FLOW) + + +def influent_at(time_day: float, weather: str) -> np.ndarray: + """Generate a deterministic BSM1 influent at a 15-minute timestamp. + + The dry profile uses the published BSM1 dynamic-load averages. Rain and + storm add dilution water over the event windows described by the IWA + report; pollutant mass flow is conserved during the dilution event. + """ + + phase = 2.0 * math.pi * (time_day % 1.0) + weekend = 0.82 if int(time_day) % 7 in (5, 6) else 1.0 + flow = 18446.0 * weekend * (1.0 + 0.31 * math.sin(phase - 1.0) + 0.08 * math.sin(2.0 * phase - 0.4)) + flow = max(9000.0, flow) + concentration_scale = 1.0 + 0.12 * math.sin(phase + 0.55) + 0.04 * math.sin(2.0 * phase) + components = np.array( + [30.0, 69.50, 51.20, 202.32, 28.17, 0.0, 0.0, 0.0, 0.0, 31.56, 6.95, 10.59, 7.0], + dtype=float, + ) + variable = [SS, XI, XS, XBH, SNH, SND, XND] + components[variable] *= concentration_scale + + extra_water = 0.0 + if weather == "rain" and 8.35 <= time_day <= 10.44: + edge = min((time_day - 8.35) / 0.08, (10.44 - time_day) / 0.08, 1.0) + extra_water = 20000.0 * max(0.0, edge) + elif weather == "storm": + extra_water = 42000.0 * math.exp(-0.5 * ((time_day - 8.87) / 0.075) ** 2) + extra_water += 45000.0 * math.exp(-0.5 * ((time_day - 11.18) / 0.16) ** 2) + if extra_water > 0.0: + dilution = flow / (flow + extra_water) + # Added water dilutes every soluble and particulate concentration, + # including alkalinity. The corresponding component mass flow is + # conserved while the hydraulic load increases. + components[:13] *= dilution + flow += extra_water + return _stream(components, flow) + + +def advanced_quantities(effluent: np.ndarray) -> dict[str, float]: + p = ASM1_PARAMETERS + kjeldahl = ( + effluent[SNH] + effluent[SND] + effluent[XND] + + p[17] * (effluent[XBH] + effluent[XBA]) + + p[18] * (effluent[XP] + effluent[XI]) + ) + cod = float(np.sum(effluent[[SI, SS, XI, XS, XBH, XBA, XP]])) + bod5 = 0.25 * (effluent[SS] + effluent[XS] + (1.0 - p[16]) * (effluent[XBH] + effluent[XBA])) + return { + "total_nitrogen": float(kjeldahl + effluent[SNO]), + "cod": cod, + "bod5": float(bod5), + "tss": float(effluent[TSS]), + } + + +def effluent_quality_index(effluent: np.ndarray) -> float: + advanced = advanced_quantities(effluent) + kjeldahl = advanced["total_nitrogen"] - effluent[SNO] + weighted = ( + 2.0 * effluent[TSS] + + advanced["cod"] + + 30.0 * kjeldahl + + 10.0 * effluent[SNO] + + 2.0 * advanced["bod5"] + ) + return float(weighted * effluent[Q] / 1000.0) + + +class BSM1Plant: + """Stateful five-reactor BSM1 plant with a ten-layer settler.""" + + def __init__(self) -> None: + self.reactors = _REACTOR_INITIAL.copy() + soluble_initial = np.array([30.0, 0.857137129, 2.0, 13.711131912, 0.866431079, 0.690377137, 3.828235655]) + self.settler = np.concatenate([np.full(SETTLER_LAYERS, value) for value in soluble_initial] + [_SETTLER_TSS_INITIAL.copy()]) + return_components = np.array([30.0, 0.857137129, 2247.12680338, 92.61753044, 5007.32236075, 299.51772354, 885.01259999, 2.0, 13.711131912, 0.866431079, 0.690377137, 6.637034575, 3.828235655]) + self.return_sludge = _stream(return_components, RETURN_FLOW) + self.internal_recycle = _stream(self.reactors[-1], 55338.0) + effluent_components = np.array([30.0, 0.857137129, 4.390198009, 0.180946308, 9.782775332, 0.585165960, 1.729043750, 2.0, 13.711131912, 0.866431079, 0.690377137, 0.012966734, 3.828235655]) + self.effluent = _stream(effluent_components, 18061.0) + + def step(self, influent: np.ndarray, action: Mapping[str, float], dt_days: float) -> np.ndarray: + qintr = float(action["internal_recycle_m3_per_day"]) + self.internal_recycle[Q] = qintr + inlet = combine_streams(influent, self.return_sludge, self.internal_recycle) + klas = [0.0, 0.0, float(action["kla3_per_day"]), float(action["kla4_per_day"]), float(action["kla5_per_day"])] + stream = inlet + for index, (kla, volume) in enumerate(zip(klas, VOLUMES)): + self.reactors[index] = _integrate_reactor(self.reactors[index], stream, kla, float(volume), dt_days) + stream = _stream(self.reactors[index], inlet[Q], influent[TEMP]) + + settler_feed = _stream(self.reactors[-1], max(stream[Q] - qintr, 1.0), influent[TEMP]) + self.internal_recycle = _stream(self.reactors[-1], qintr, influent[TEMP]) + self.settler = odeint( + lambda y, _: _settler_rhs(y, settler_feed), + self.settler, + np.array([0.0, dt_days]), + rtol=1e-4, + atol=1e-6, + mxstep=1000, + )[-1] + if not np.all(np.isfinite(self.settler)): + raise FloatingPointError("settler integration produced non-finite state") + self.settler = np.maximum(self.settler, 0.0) + self.effluent, self.return_sludge = _settler_outputs(self.settler, settler_feed) + return self.effluent.copy() + + +def action_energy(action: Mapping[str, float]) -> dict[str, float]: + klas = np.array([0.0, 0.0, action["kla3_per_day"], action["kla4_per_day"], action["kla5_per_day"]], dtype=float) + aeration = float(np.sum(OXYGEN_SATURATION * VOLUMES * klas) / (1.8 * 1000.0)) + pumping = float(0.004 * action["internal_recycle_m3_per_day"] + 0.008 * RETURN_FLOW + 0.05 * WASTE_FLOW) + mixing = float(24.0 * 0.005 * np.sum(VOLUMES[klas < 20.0])) + return {"aeration": aeration, "pumping": pumping, "mixing": mixing} diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/evaluator.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/evaluator.py new file mode 100644 index 00000000..c2cdd3ac --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/evaluator.py @@ -0,0 +1,309 @@ +"""Evaluator for the BSM1AerationControl benchmark.""" + +from __future__ import annotations + +import argparse +import json +import math +from pathlib import Path +from typing import Any, Callable, Mapping + +import numpy as np + +try: + from .bsm1_model import ( + BSM1Plant, SNH, SNO, SO, action_energy, advanced_quantities, + effluent_quality_index, influent_at, load_config, + ) + from .policy_runtime import PolicyRuntime +except ImportError: + from bsm1_model import ( + BSM1Plant, SNH, SNO, SO, action_energy, advanced_quantities, + effluent_quality_index, influent_at, load_config, + ) + from policy_runtime import PolicyRuntime + + +ROOT = Path(__file__).resolve().parents[1] +DEFAULT_CANDIDATE = ROOT / "scripts" / "init.py" +ACTION_KEYS = ( + "kla3_per_day", "kla4_per_day", "kla5_per_day", "internal_recycle_m3_per_day", +) + + +def validate_action(value: Any, previous: Mapping[str, float], config: Mapping[str, Any]) -> dict[str, float]: + if not isinstance(value, dict): + raise TypeError("control must return a JSON object") + if set(value) != set(ACTION_KEYS): + missing = sorted(set(ACTION_KEYS) - set(value)) + extra = sorted(set(value) - set(ACTION_KEYS)) + raise ValueError(f"action keys mismatch; missing={missing}, extra={extra}") + action: dict[str, float] = {} + for key in ACTION_KEYS: + raw = value[key] + if isinstance(raw, bool) or not isinstance(raw, (int, float)): + raise TypeError(f"action[{key!r}] must be a finite number") + number = float(raw) + if not math.isfinite(number): + raise ValueError(f"action[{key!r}] must be finite") + action[key] = number + + limits = config["action"] + for key in ACTION_KEYS[:3]: + if not limits["kla_min_per_day"] <= action[key] <= limits["kla_max_per_day"]: + raise ValueError(f"action[{key!r}] is outside the KLa range") + if abs(action[key] - previous[key]) > limits["kla_max_change_per_step"] + 1e-9: + raise ValueError(f"action[{key!r}] exceeds the per-step KLa slew limit") + recycle_key = "internal_recycle_m3_per_day" + if not limits["qintr_min_m3_per_day"] <= action[recycle_key] <= limits["qintr_max_m3_per_day"]: + raise ValueError("internal recycle flow is outside its range") + if abs(action[recycle_key] - previous[recycle_key]) > limits["qintr_max_change_per_step"] + 1e-9: + raise ValueError("internal recycle flow exceeds its per-step slew limit") + return action + + +def _observation( + plant: BSM1Plant, + scenario: Mapping[str, str], + time_day: float, + influent: np.ndarray, + previous: Mapping[str, float], + noise: Mapping[str, float], + step_minutes: float, +) -> dict[str, Any]: + advanced = advanced_quantities(plant.effluent) + return { + "scenario_id": scenario["scenario_id"], + "weather": scenario["weather"], + "time_day": float(time_day), + "step_minutes": float(step_minutes), + "influent_flow_m3_per_day": float(influent[14]), + "influent_ammonium_gN_per_m3": float(influent[SNH]), + "reactor3_do_gO2_per_m3": max(0.0, float(plant.reactors[2, SO]) + noise["do3"]), + "reactor4_do_gO2_per_m3": max(0.0, float(plant.reactors[3, SO]) + noise["do4"]), + "reactor5_do_gO2_per_m3": max(0.0, float(plant.reactors[4, SO]) + noise["do5"]), + "reactor2_nitrate_gN_per_m3": max(0.0, float(plant.reactors[1, SNO]) + noise["nitrate"]), + "effluent_ammonium_gN_per_m3": max(0.0, float(plant.effluent[SNH]) + noise["ammonium"]), + "effluent_total_nitrogen_gN_per_m3": float(advanced["total_nitrogen"]), + "previous_action": dict(previous), + } + + +def _score(metrics: Mapping[str, float], config: Mapping[str, Any]) -> float: + scoring = config["scoring"] + quality = math.exp( + -metrics["eqi_kg_pollution_units_per_day"] / scoring["eqi_scale"] + ) + aeration = math.exp( + -metrics["aeration_energy_kwh_per_day"] / scoring["aeration_energy_scale"] + ) + pumping = math.exp( + -metrics["pumping_energy_kwh_per_day"] / scoring["pumping_energy_scale"] + ) + mixing = math.exp( + -metrics["mixing_energy_kwh_per_day"] / scoring["mixing_energy_scale"] + ) + smoothness = math.exp(-scoring["switching_exponent"] * metrics["switching_index"]) + weights = scoring["objective_weights"] + base = 100.0 * ( + weights["quality"] * quality + + weights["aeration"] * aeration + + weights["pumping"] * pumping + + weights["mixing"] * mixing + + weights["smoothness"] * smoothness + ) + compliance = math.exp( + -scoring["violation_exponent"] * metrics["violation_index"] + ) + return float(np.clip(base * compliance, 0.0, 100.0)) + + +def simulate_scenario( + scenario: Mapping[str, str], + action_provider: Callable[[dict[str, Any]], Any], + reset_provider: Callable[[dict[str, Any]], None], + *, + max_steps: int | None = None, +) -> dict[str, Any]: + config = load_config() + simulation = config["simulation"] + step_minutes = float(simulation["step_minutes"]) + dt_days = step_minutes / (24.0 * 60.0) + total_steps = int(round(float(simulation["days"]) / dt_days)) + if max_steps is not None: + total_steps = min(total_steps, int(max_steps)) + evaluation_start = float(simulation["evaluation_start_day"]) + if max_steps is not None and total_steps * dt_days <= evaluation_start: + evaluation_start = 0.0 + + descriptor = { + "scenario_id": scenario["scenario_id"], "weather": scenario["weather"], + "duration_days": total_steps * dt_days, "step_minutes": step_minutes, + } + reset_provider(descriptor) + plant = BSM1Plant() + previous = {key: float(config["reference_action"][key]) for key in ACTION_KEYS} + seed_offset = {"dry": 0, "rain": 1000, "storm": 2000}[scenario["weather"]] + rng = np.random.default_rng(int(simulation["sensor_noise_seed"]) + seed_offset) + std = simulation["sensor_noise_std"] + + eqi_values: list[float] = [] + aeration_values: list[float] = [] + pumping_values: list[float] = [] + mixing_values: list[float] = [] + switching_values: list[float] = [] + violation_values: list[float] = [] + violation_excess: dict[str, list[float]] = { + name: [] for name in ("ammonium", "total_nitrogen", "cod", "tss", "bod5") + } + maximum_quantities = {name: 0.0 for name in violation_excess} + effluent_rows: list[dict[str, float]] = [] + + for index in range(total_steps): + time_day = index * dt_days + influent = influent_at(time_day, scenario["weather"]) + noise = { + "do3": rng.normal(0.0, std["dissolved_oxygen"]), + "do4": rng.normal(0.0, std["dissolved_oxygen"]), + "do5": rng.normal(0.0, std["dissolved_oxygen"]), + "nitrate": rng.normal(0.0, std["nitrate"]), + "ammonium": rng.normal(0.0, std["ammonium"]), + } + observation = _observation(plant, scenario, time_day, influent, previous, noise, step_minutes) + action = validate_action(action_provider(observation), previous, config) + effluent = plant.step(influent, action, dt_days) + if time_day >= evaluation_start: + advanced = advanced_quantities(effluent) + energy = action_energy(action) + limits = config["limits"] + quantities = { + "ammonium": float(effluent[SNH]), + "total_nitrogen": advanced["total_nitrogen"], + "cod": advanced["cod"], "tss": advanced["tss"], "bod5": advanced["bod5"], + } + limit_keys = { + "ammonium": "ammonium_gN_per_m3", + "total_nitrogen": "total_nitrogen_gN_per_m3", + "cod": "cod_gCOD_per_m3", + "tss": "tss_gSS_per_m3", + "bod5": "bod5_gBOD_per_m3", + } + normalized_excess = [ + max(0.0, quantities[name] / limits[limit_key] - 1.0) + for name, limit_key in limit_keys.items() + ] + for name, excess in zip(limit_keys, normalized_excess): + violation_excess[name].append(float(excess)) + maximum_quantities[name] = max( + maximum_quantities[name], float(quantities[name]) + ) + switching = sum( + abs(action[key] - previous[key]) / (360.0 if key.startswith("kla") else 92230.0) + for key in ACTION_KEYS + ) + eqi_values.append(effluent_quality_index(effluent)) + aeration_values.append(energy["aeration"]) + pumping_values.append(energy["pumping"]) + mixing_values.append(energy["mixing"]) + switching_values.append(float(switching)) + violation_values.append(float(np.mean(normalized_excess))) + samples_per_day = max(1, int(round(24.0 * 60.0 / step_minutes))) + if index % samples_per_day == 0 or index == total_steps - 1: + effluent_rows.append({"time_day": time_day, **quantities}) + previous = action + + if not eqi_values: + raise ValueError("scenario has no evaluation samples") + metrics = { + "eqi_kg_pollution_units_per_day": float(np.mean(eqi_values)), + "aeration_energy_kwh_per_day": float(np.mean(aeration_values)), + "pumping_energy_kwh_per_day": float(np.mean(pumping_values)), + "mixing_energy_kwh_per_day": float(np.mean(mixing_values)), + "switching_index": float(np.mean(switching_values)), + "violation_index": float(np.mean(violation_values)), + "max_ammonium_gN_per_m3": maximum_quantities["ammonium"], + "max_total_nitrogen_gN_per_m3": maximum_quantities["total_nitrogen"], + "max_cod_gCOD_per_m3": maximum_quantities["cod"], + "max_tss_gSS_per_m3": maximum_quantities["tss"], + "max_bod5_gBOD_per_m3": maximum_quantities["bod5"], + } + for name, values in violation_excess.items(): + metrics[f"{name}_violation_fraction"] = float(np.mean(np.asarray(values) > 0.0)) + metrics[f"{name}_mean_normalized_excess"] = float(np.mean(values)) + return { + "scenario_id": scenario["scenario_id"], "weather": scenario["weather"], + "score": _score(metrics, config), "metrics": metrics, "effluent_samples": effluent_rows, + } + + +def evaluate(candidate_path: Path, *, max_steps: int | None = None) -> dict[str, Any]: + config = load_config() + rows: list[dict[str, Any]] = [] + try: + with PolicyRuntime(candidate_path) as runtime: + for scenario in config["scenarios"]: + rows.append(simulate_scenario( + scenario, runtime.control, runtime.reset, max_steps=max_steps, + )) + except Exception as exc: + return { + "combined_score": 0.0, "diagnostic_score": 0.0, "valid": 0.0, + "completed_scenarios": float(len(rows)), "error": f"{type(exc).__name__}: {exc}", + "rows": rows, + } + scores = np.asarray([row["score"] for row in rows], dtype=float) + scoring = config["scoring"] + diagnostic = float( + scoring["mean_weight"] * np.mean(scores) + + scoring["worst_case_weight"] * np.min(scores) + ) + return { + "combined_score": diagnostic, "diagnostic_score": diagnostic, "valid": 1.0, + "completed_scenarios": float(len(rows)), "rows": rows, + } + + +def _write_json(path: str | None, payload: Mapping[str, Any]) -> None: + if not path: + return + output = Path(path).expanduser().resolve() + output.parent.mkdir(parents=True, exist_ok=True) + output.write_text(json.dumps(payload, indent=2, ensure_ascii=False), encoding="utf-8") + + +def main() -> None: + parser = argparse.ArgumentParser(description="Evaluate a BSM1 aeration controller") + parser.add_argument("candidate", nargs="?", default=str(DEFAULT_CANDIDATE)) + parser.add_argument("--metrics-out", default=None) + parser.add_argument("--artifacts-out", default=None) + args = parser.parse_args() + candidate = Path(args.candidate).expanduser().resolve() + result = evaluate(candidate) + print("=== BSM1 Aeration Control ===") + for row in result["rows"]: + metrics = row["metrics"] + print( + f"scenario={row['scenario_id']} score={row['score']:.3f} " + f"eqi={metrics['eqi_kg_pollution_units_per_day']:.2f} " + f"ae={metrics['aeration_energy_kwh_per_day']:.2f} " + f"violations={metrics['violation_index']:.6f}" + ) + if "error" in result: + print(f"error: {result['error']}") + print("---") + print(f"completed_scenarios: {result['completed_scenarios']:.0f}/{len(load_config()['scenarios'])}") + print(f"diagnostic_score: {result['diagnostic_score']:.4f}") + print(f"combined_score: {result['combined_score']:.4f}") + metrics_out = { + "combined_score": float(result["combined_score"]), + "diagnostic_score": float(result["diagnostic_score"]), + "valid": float(result["valid"]), + "completed_scenarios": float(result["completed_scenarios"]), + "num_scenarios": float(len(load_config()["scenarios"])), + } + _write_json(args.metrics_out, metrics_out) + _write_json(args.artifacts_out, {"candidate_path": candidate.name, "rows": result["rows"]}) + + +if __name__ == "__main__": + main() diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_runtime.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_runtime.py new file mode 100644 index 00000000..9227f7b8 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_runtime.py @@ -0,0 +1,251 @@ +"""Bounded parent-side runtime for BSM1 candidate controllers.""" + +from __future__ import annotations + +import json +import os +import queue +import shutil +import signal +import subprocess +import sys +import tempfile +import threading +import time +from collections import deque +from pathlib import Path +from typing import Any + + +MAX_RESPONSE_BYTES = 16 * 1024 +_EOF = object() + + +class PolicyRuntimeError(RuntimeError): + pass + + +class PolicyTimeoutError(PolicyRuntimeError): + pass + + +class PolicyProtocolError(PolicyRuntimeError): + pass + + +class PolicyCandidateError(PolicyRuntimeError): + pass + + +def _clean_environment() -> dict[str, str]: + allowed = { + "COMSPEC", "LANG", "LC_ALL", "NUMBER_OF_PROCESSORS", "PATH", "PATHEXT", + "SYSTEMDRIVE", "SYSTEMROOT", "TEMP", "TMP", "WINDIR", + } + env = {key: value for key, value in os.environ.items() if key.upper() in allowed} + env["PYTHONUTF8"] = "1" + env["PYTHONIOENCODING"] = "utf-8" + return env + + +class PolicyRuntime: + def __init__( + self, + candidate_path: str | Path, + *, + startup_timeout_s: float = 3.0, + call_timeout_s: float = 0.10, + candidate_time_budget_s: float = 12.0, + ) -> None: + if min(startup_timeout_s, call_timeout_s, candidate_time_budget_s) <= 0: + raise ValueError("policy timeouts must be positive") + source = Path(candidate_path).expanduser().resolve() + if not source.is_file(): + raise FileNotFoundError(source) + self.call_timeout_s = float(call_timeout_s) + self._candidate_time_remaining = float(candidate_time_budget_s) + self._next_id = 1 + self._closed = False + self._responses: queue.Queue[Any] = queue.Queue() + self._stderr_chunks: deque[bytes] = deque() + self._stderr_size = 0 + self._tempdir = Path(tempfile.mkdtemp(prefix="bsm1_policy_")) + self._process: subprocess.Popen[bytes] | None = None + try: + copied_candidate = self._tempdir / "candidate.py" + shutil.copy2(source, copied_candidate) + worker = Path(__file__).with_name("policy_worker.py").resolve() + kwargs: dict[str, Any] = { + "cwd": str(self._tempdir), "env": _clean_environment(), "stdin": subprocess.PIPE, + "stdout": subprocess.PIPE, "stderr": subprocess.PIPE, "bufsize": 0, + } + if os.name == "nt": + kwargs["creationflags"] = subprocess.CREATE_NEW_PROCESS_GROUP + else: + kwargs["start_new_session"] = True + self._process = subprocess.Popen( + [sys.executable, "-I", "-u", str(worker), str(copied_candidate)], **kwargs + ) + assert self._process.stdout is not None and self._process.stderr is not None + threading.Thread(target=self._read_stdout, args=(self._process.stdout,), daemon=True).start() + threading.Thread(target=self._drain_stderr, args=(self._process.stderr,), daemon=True).start() + ready = self._receive(startup_timeout_s, "candidate import", charge_budget=False) + if not isinstance(ready, dict) or not ready.get("ok") or not ready.get("ready"): + self._raise_response_error(ready, "candidate import") + except BaseException: + self.close(force=True) + raise + + @property + def stderr_tail(self) -> str: + return b"".join(self._stderr_chunks).decode("utf-8", errors="replace") + + def _read_stdout(self, stream: Any) -> None: + try: + while True: + line = stream.readline(MAX_RESPONSE_BYTES + 1) + if not line: + break + self._responses.put(line) + if len(line) > MAX_RESPONSE_BYTES or not line.endswith(b"\n"): + break + finally: + self._responses.put(_EOF) + + def _drain_stderr(self, stream: Any) -> None: + while True: + chunk = stream.read(4096) + if not chunk: + return + self._stderr_chunks.append(chunk) + self._stderr_size += len(chunk) + while self._stderr_size > 16 * 1024 and self._stderr_chunks: + self._stderr_size -= len(self._stderr_chunks.popleft()) + + def _receive(self, timeout_s: float, operation: str, *, charge_budget: bool = True) -> dict[str, Any]: + if charge_budget and self._candidate_time_remaining <= 0.0: + self.close(force=True) + raise PolicyTimeoutError("candidate cumulative runtime budget exceeded") + timeout = min(timeout_s, self._candidate_time_remaining) if charge_budget else timeout_s + started = time.monotonic() + try: + item = self._responses.get(timeout=timeout) + except queue.Empty as exc: + self.close(force=True) + raise PolicyTimeoutError(f"{operation} timed out after {timeout:.3f}s") from exc + finally: + if charge_budget: + self._candidate_time_remaining -= time.monotonic() - started + if item is _EOF: + code = None if self._process is None else self._process.poll() + detail = self.stderr_tail[-1000:] + raise PolicyCandidateError( + f"candidate worker exited unexpectedly (code={code})" + (f": {detail}" if detail else "") + ) + if not isinstance(item, bytes) or len(item) > MAX_RESPONSE_BYTES or not item.endswith(b"\n"): + self.close(force=True) + raise PolicyProtocolError("candidate response exceeds 16 KiB or lacks newline") + try: + response = json.loads(item.decode("utf-8")) + except (UnicodeDecodeError, json.JSONDecodeError) as exc: + self.close(force=True) + raise PolicyProtocolError("candidate returned invalid UTF-8 JSON") from exc + if not isinstance(response, dict): + raise PolicyProtocolError("candidate response must be a JSON object") + return response + + @staticmethod + def _raise_response_error(response: Any, operation: str) -> None: + if isinstance(response, dict) and isinstance(response.get("error"), dict): + error = response["error"] + raise PolicyCandidateError( + f"{operation} failed: {error.get('type', 'CandidateError')}: {error.get('message', '')}" + ) + raise PolicyProtocolError(f"malformed response during {operation}") + + def _rpc(self, operation: str, **payload: Any) -> dict[str, Any]: + if self._closed or self._process is None or self._process.poll() is not None: + raise PolicyCandidateError("candidate worker is not running") + request_id = self._next_id + self._next_id += 1 + encoded = (json.dumps( + {"id": request_id, "op": operation, **payload}, ensure_ascii=True, + allow_nan=False, separators=(",", ":"), + ) + "\n").encode("utf-8") + try: + assert self._process.stdin is not None + self._process.stdin.write(encoded) + self._process.stdin.flush() + except (BrokenPipeError, OSError) as exc: + raise PolicyCandidateError("candidate worker closed its input") from exc + response = self._receive(self.call_timeout_s, operation) + if response.get("id") != request_id: + self.close(force=True) + raise PolicyProtocolError("candidate response id does not match request") + if not response.get("ok"): + self._raise_response_error(response, operation) + return response + + def reset(self, scenario: dict[str, Any]) -> None: + self._rpc("reset", scenario=scenario) + + def control(self, observation: dict[str, Any]) -> Any: + return self._rpc("control", observation=observation).get("action") + + def close(self, *, force: bool = False) -> None: + if self._closed: + return + process = self._process + if process is not None and process.poll() is None and not force: + try: + self._rpc("shutdown") + process.wait(timeout=0.25) + except Exception: + force = True + if process is not None and process.poll() is None: + self._kill_process_tree(process) + self._closed = True + if process is not None: + for stream in (process.stdin, process.stdout, process.stderr): + try: + if stream is not None: + stream.close() + except OSError: + pass + shutil.rmtree(self._tempdir, ignore_errors=True) + + @staticmethod + def _kill_process_tree(process: subprocess.Popen[bytes]) -> None: + if process.poll() is not None: + return + if os.name == "nt": + try: + subprocess.run( + ["taskkill", "/PID", str(process.pid), "/T", "/F"], + stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, + timeout=2.0, check=False, + ) + except Exception: + process.kill() + else: + try: + os.killpg(process.pid, signal.SIGKILL) + except (ProcessLookupError, PermissionError): + process.kill() + try: + process.wait(timeout=2.0) + except subprocess.TimeoutExpired: + process.kill() + process.wait(timeout=2.0) + + def __enter__(self) -> "PolicyRuntime": + return self + + def __exit__(self, exc_type: Any, exc: Any, traceback: Any) -> None: + self.close() + + def __del__(self) -> None: + try: + self.close(force=True) + except Exception: + pass diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py new file mode 100644 index 00000000..2abaec4f --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py @@ -0,0 +1,98 @@ +"""JSON-lines worker for an untrusted BSM1 controller module. + +This provides process and protocol isolation, not an operating-system sandbox. +""" + +from __future__ import annotations + +import importlib.util +import json +import os +import sys +from pathlib import Path +from typing import Any + + +_json_loads = json.loads +_json_dumps = json.dumps + + +def _protocol_streams() -> tuple[Any, Any]: + protocol_in = os.fdopen(os.dup(0), "r", encoding="utf-8", newline="\n") + protocol_out = os.fdopen(os.dup(1), "w", encoding="utf-8", newline="\n", buffering=1) + devnull = os.open(os.devnull, os.O_WRONLY) + try: + os.dup2(devnull, 1) + finally: + os.close(devnull) + return protocol_in, protocol_out + + +def _load_candidate(path: Path) -> Any: + spec = importlib.util.spec_from_file_location("bsm1_candidate", path) + if spec is None or spec.loader is None: + raise ImportError(f"cannot load candidate from {path}") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + if not callable(getattr(module, "control", None)): + raise AttributeError("candidate must define callable control(observation)") + if not callable(getattr(module, "reset_controller", None)): + raise AttributeError("candidate must define callable reset_controller(scenario)") + return module + + +def _error(request_id: Any, exc: BaseException) -> dict[str, Any]: + try: + message = str(exc)[:500] + except Exception: + message = "failed to format candidate exception" + return {"id": request_id, "ok": False, "error": {"type": type(exc).__name__, "message": message}} + + +def _send(stream: Any, payload: dict[str, Any]) -> None: + stream.write(_json_dumps(payload, ensure_ascii=True, allow_nan=False, separators=(",", ":")) + "\n") + stream.flush() + + +def main() -> int: + protocol_in, protocol_out = _protocol_streams() + if len(sys.argv) != 2: + _send(protocol_out, _error(None, ValueError("expected candidate path"))) + return 2 + try: + candidate = _load_candidate(Path(sys.argv[1]).resolve()) + except BaseException as exc: + _send(protocol_out, _error(None, exc)) + return 1 + _send(protocol_out, {"id": None, "ok": True, "ready": True}) + + for line in protocol_in: + request_id: Any = None + try: + request = _json_loads(line) + if not isinstance(request, dict): + raise TypeError("request must be a JSON object") + request_id = request.get("id") + operation = request.get("op") + if operation == "reset": + candidate.reset_controller(request.get("scenario")) + response = {"id": request_id, "ok": True} + elif operation == "control": + action = candidate.control(request.get("observation")) + response = {"id": request_id, "ok": True, "action": action} + elif operation == "shutdown": + _send(protocol_out, {"id": request_id, "ok": True}) + return 0 + else: + raise ValueError(f"unknown operation: {operation!r}") + _send(protocol_out, response) + except BaseException as exc: + try: + _send(protocol_out, _error(request_id, exc)) + except BaseException: + return 1 + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/requirements.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/requirements.txt new file mode 100644 index 00000000..0ffddefd --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/requirements.txt @@ -0,0 +1,2 @@ +numpy>=1.24 +scipy>=1.10 diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py new file mode 100644 index 00000000..864f2bba --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py @@ -0,0 +1,50 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +try: + from verification.bsm1_model import ( + Q, + BSM1Plant, + action_energy, + influent_at, + ) +except ModuleNotFoundError: + from bsm1_model import Q, BSM1Plant, action_energy, influent_at + + +BASELINE_ACTION = { + "kla3_per_day": 240.0, + "kla4_per_day": 240.0, + "kla5_per_day": 84.0, + "internal_recycle_m3_per_day": 55338.0, +} + + +class BSM1ModelTests(unittest.TestCase): + def test_reference_energy_matches_published_bsm1_operating_point(self) -> None: + energy = action_energy(BASELINE_ACTION) + self.assertAlmostEqual(energy["aeration"], 3341.3866666666663, places=9) + self.assertAlmostEqual(energy["pumping"], 388.17, places=9) + self.assertAlmostEqual(energy["mixing"], 240.0, places=9) + + def test_rain_adds_water_without_changing_component_mass_flow(self) -> None: + time_day = 9.0 + dry = influent_at(time_day, "dry") + rain = influent_at(time_day, "rain") + self.assertAlmostEqual(rain[Q] - dry[Q], 20000.0, places=9) + np.testing.assert_allclose(rain[:13] * rain[Q], dry[:13] * dry[Q], rtol=1e-12) + + def test_recycle_change_keeps_the_settler_hydraulically_valid(self) -> None: + plant = BSM1Plant() + action = dict(BASELINE_ACTION, internal_recycle_m3_per_day=30000.0) + effluent = plant.step(influent_at(0.0, "dry"), action, 15.0 / 1440.0) + self.assertTrue(np.all(np.isfinite(effluent))) + self.assertGreater(effluent[Q], 0.0) + self.assertAlmostEqual(effluent[Q], influent_at(0.0, "dry")[Q] - 385.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_evaluator.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_evaluator.py new file mode 100644 index 00000000..2f745453 --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_evaluator.py @@ -0,0 +1,81 @@ +from __future__ import annotations + +import tempfile +import unittest +from pathlib import Path +from typing import Any + +try: + from verification.bsm1_model import load_config + from verification.evaluator import evaluate, simulate_scenario, validate_action +except ModuleNotFoundError: + from bsm1_model import load_config + from evaluator import evaluate, simulate_scenario, validate_action + + +BASELINE_ACTION = { + "kla3_per_day": 240.0, + "kla4_per_day": 240.0, + "kla5_per_day": 84.0, + "internal_recycle_m3_per_day": 55338.0, +} + + +class EvaluatorTests(unittest.TestCase): + def _candidate(self, source: str) -> Path: + directory = tempfile.TemporaryDirectory(prefix="bsm1_candidate_test_") + self.addCleanup(directory.cleanup) + path = Path(directory.name) / "candidate.py" + path.write_text(source, encoding="utf-8") + return path + + def test_action_schema_type_range_and_slew_are_enforced(self) -> None: + config = load_config() + self.assertEqual(validate_action(dict(BASELINE_ACTION), BASELINE_ACTION, config), BASELINE_ACTION) + invalid_actions: list[dict[str, Any]] = [ + {key: value for key, value in BASELINE_ACTION.items() if key != "kla5_per_day"}, + dict(BASELINE_ACTION, kla3_per_day=True), + dict(BASELINE_ACTION, kla3_per_day=361.0), + dict(BASELINE_ACTION, internal_recycle_m3_per_day=1000.0), + ] + for action in invalid_actions: + with self.subTest(action=action), self.assertRaises((TypeError, ValueError)): + validate_action(action, BASELINE_ACTION, config) + + def test_lower_recycle_exposes_measurable_optimization_headroom(self) -> None: + scenario = {"scenario_id": "dry", "weather": "dry"} + lower_recycle = dict(BASELINE_ACTION, internal_recycle_m3_per_day=30000.0) + reset = lambda _: None + baseline = simulate_scenario( + scenario, lambda _: BASELINE_ACTION, reset, max_steps=96 + ) + improved = simulate_scenario( + scenario, lambda _: lower_recycle, reset, max_steps=96 + ) + self.assertGreater(improved["score"], baseline["score"] + 0.5) + self.assertLess( + improved["metrics"]["pumping_energy_kwh_per_day"], + baseline["metrics"]["pumping_energy_kwh_per_day"], + ) + + def test_short_baseline_evaluation_is_valid_and_deterministic(self) -> None: + source = ( + "def reset_controller(scenario): pass\n" + f"def control(observation): return {BASELINE_ACTION!r}\n" + ) + candidate = self._candidate(source) + first = evaluate(candidate, max_steps=2) + second = evaluate(candidate, max_steps=2) + self.assertEqual(first, second) + self.assertEqual(first["valid"], 1.0) + self.assertEqual(first["completed_scenarios"], 3.0) + + def test_malformed_candidate_returns_zero_instead_of_crashing(self) -> None: + result = evaluate(self._candidate("def broken(:\n"), max_steps=2) + self.assertEqual(result["valid"], 0.0) + self.assertEqual(result["combined_score"], 0.0) + self.assertIn("candidate import failed", result["error"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py new file mode 100644 index 00000000..a3661a8e --- /dev/null +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py @@ -0,0 +1,46 @@ +from __future__ import annotations + +import tempfile +import unittest +from pathlib import Path + +try: + from verification.policy_runtime import PolicyRuntime, PolicyTimeoutError +except ModuleNotFoundError: + from policy_runtime import PolicyRuntime, PolicyTimeoutError + + +class PolicyRuntimeTests(unittest.TestCase): + def _candidate(self, source: str) -> Path: + directory = tempfile.TemporaryDirectory(prefix="bsm1_runtime_test_") + self.addCleanup(directory.cleanup) + path = Path(directory.name) / "candidate.py" + path.write_text(source, encoding="utf-8") + return path + + def test_stdout_is_isolated_from_the_json_protocol(self) -> None: + candidate = self._candidate( + "print('during import')\n" + "def reset_controller(scenario): print('during reset')\n" + "def control(observation):\n" + " print('during control')\n" + " return {'ok': 1}\n" + ) + with PolicyRuntime(candidate) as runtime: + runtime.reset({"scenario_id": "dry"}) + self.assertEqual(runtime.control({"step": 0}), {"ok": 1}) + + def test_slow_control_is_terminated(self) -> None: + candidate = self._candidate( + "import time\n" + "def reset_controller(scenario): pass\n" + "def control(observation): time.sleep(0.2); return {}\n" + ) + with PolicyRuntime(candidate, call_timeout_s=0.03) as runtime: + runtime.reset({}) + with self.assertRaises(PolicyTimeoutError): + runtime.control({}) + + +if __name__ == "__main__": + unittest.main() diff --git a/benchmarks/WastewaterTreatment/README.md b/benchmarks/WastewaterTreatment/README.md new file mode 100644 index 00000000..d500038a --- /dev/null +++ b/benchmarks/WastewaterTreatment/README.md @@ -0,0 +1,11 @@ +# Wastewater Treatment + +This domain collects control and optimization tasks for biological wastewater-treatment +plants. The benchmarks emphasize process dynamics, effluent compliance, energy use, and +robust operation under changing hydraulic loads. + +## Tasks + +- `BSM1AerationControl` + - Unified benchmark: `task=unified task.benchmark=WastewaterTreatment/BSM1AerationControl` + - Description: optimize zone aeration and internal recycle in an IWA BSM1-derived activated-sludge plant across dry, rain, and storm operation. diff --git a/benchmarks/WastewaterTreatment/README_zh-CN.md b/benchmarks/WastewaterTreatment/README_zh-CN.md new file mode 100644 index 00000000..a6880660 --- /dev/null +++ b/benchmarks/WastewaterTreatment/README_zh-CN.md @@ -0,0 +1,9 @@ +# 污水处理 + +本领域收集生物污水处理厂的控制与优化任务，重点关注过程动态、出水达标、能源消耗以及水力负荷变化下的稳健运行。 + +## 任务列表 + +- `BSM1AerationControl` + - 统一任务：`task=unified task.benchmark=WastewaterTreatment/BSM1AerationControl` + - 简介：在基于 IWA BSM1 的活性污泥厂中，针对旱天、降雨和暴雨工况优化分区曝气与内回流。 diff --git a/frontier_eval/conf/task/bsm1_aeration_control.yaml b/frontier_eval/conf/task/bsm1_aeration_control.yaml new file mode 100644 index 00000000..0751974e --- /dev/null +++ b/frontier_eval/conf/task/bsm1_aeration_control.yaml @@ -0,0 +1,6 @@ +name: unified +benchmark: WastewaterTreatment/BSM1AerationControl + +runtime: + env: + PYTHONUTF8: "1" From 2d5d01676091a1798fa144ed9da1a00e96dcec92 Mon Sep 17 00:00:00 2001 From: Dan Martin <93075337+Mingqwqqaq@users.noreply.github.com> Date: Wed, 19 Aug 2026 13:26:21 +0800 Subject: [PATCH 2/2] fix(WastewaterTreatment): address BSM1 review findings --- .../baseline/result_log.txt | 4 +- .../references/design_notes.md | 7 +- .../verification/policy_worker.py | 13 ++- .../verification/test_bsm1_model.py | 101 +++++++++++++++++- .../verification/test_policy_runtime.py | 15 +++ 5 files changed, 133 insertions(+), 7 deletions(-) diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt b/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt index 954cc0da..05dc79cf 100644 --- a/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/baseline/result_log.txt @@ -2,7 +2,7 @@ Direct command: python verification/evaluator.py scripts/init.py --metrics-out metrics.json --artifacts-out artifacts.json Environment: -Windows 11, Python 3.11, CPU-only +Windows 11, Python 3.12.13, CPU-only Result: === BSM1 Aeration Control === @@ -16,7 +16,7 @@ combined_score: 47.7973 Regression: python -m unittest discover -s verification -p "test_*.py" -v -Ran 9 tests in 3.574s +Ran 12 tests in 3.266s OK Unified evaluation: diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md index 10855221..0237bd6c 100644 --- a/benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/references/design_notes.md @@ -26,8 +26,11 @@ not as an exact replay of the official influent files. - The fixed operating point yields 3341.3867 kWh/day aeration energy, matching the BSM1 reference equation and published operating point. -- Unit tests check weather-event component mass conservation, positive settler hydraulics after a - recycle-flow change, action validation, deterministic results, worker isolation, and timeouts. +- Unit tests compare ASM1 and Takacs derivatives with independent numerical oracles generated by + [`bsm2-python` commit `73caa8b`](https://github.com/fau-evt/bsm2-python/commit/73caa8bde193c0dd8a8803dfcd67b546501c950e), + and also check weather-event component mass conservation, + positive settler hydraulics after a recycle-flow change, action validation, deterministic + results, worker isolation, standard dataclass imports, and timeouts. - A 96-step dry-weather regression demonstrates optimization headroom: lowering internal recycle reduces pumping energy and improves score relative to the starter policy. - The full direct evaluation covers all three 14-day scenarios; its recorded output is in diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py index 2abaec4f..be3509b5 100644 --- a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/policy_worker.py @@ -29,11 +29,20 @@ def _protocol_streams() -> tuple[Any, Any]: def _load_candidate(path: Path) -> Any: - spec = importlib.util.spec_from_file_location("bsm1_candidate", path) + module_name = "bsm1_candidate" + spec = importlib.util.spec_from_file_location(module_name, path) if spec is None or spec.loader is None: raise ImportError(f"cannot load candidate from {path}") module = importlib.util.module_from_spec(spec) - spec.loader.exec_module(module) + # exec_module() does not register manually created modules. Standard + # library features such as dataclasses with postponed annotations resolve + # types through sys.modules while the module is executing. + sys.modules[module_name] = module + try: + spec.loader.exec_module(module) + except BaseException: + sys.modules.pop(module_name, None) + raise if not callable(getattr(module, "control", None)): raise AttributeError("candidate must define callable control(observation)") if not callable(getattr(module, "reset_controller", None)): diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py index 864f2bba..720ccc7c 100644 --- a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_bsm1_model.py @@ -8,11 +8,22 @@ from verification.bsm1_model import ( Q, BSM1Plant, + _asm1_rhs, + _settler_rhs, + _stream, action_energy, influent_at, ) except ModuleNotFoundError: - from bsm1_model import Q, BSM1Plant, action_energy, influent_at + from bsm1_model import ( + Q, + BSM1Plant, + _asm1_rhs, + _settler_rhs, + _stream, + action_energy, + influent_at, + ) BASELINE_ACTION = { @@ -24,6 +35,94 @@ class BSM1ModelTests(unittest.TestCase): + # These derivative oracles were generated independently with the + # BSD-3-Clause bsm2-python equations at commit 73caa8b. Unlike invariant + # checks, they detect stoichiometric, component-mapping, and flux-sign + # regressions in the adapted implementation. + def test_asm1_derivative_matches_bsm2_python_oracle(self) -> None: + state = np.array( + [30.0, 2.5, 1100.0, 75.0, 2500.0, 150.0, 440.0, + 1.7, 9.5, 4.2, 1.0, 4.5, 4.3] + ) + feed = _stream( + np.array( + [30.0, 20.0, 100.0, 80.0, 300.0, 10.0, 50.0, + 0.5, 5.0, 20.0, 4.0, 8.0, 6.0] + ), + 75000.0, + 15.0, + ) + expected = np.array( + [ + 0.0, + -237.99652553781368, + -56264.066016504126, + -708.8242043304347, + -122581.47155209856, + -7835.430780772116, + -21882.385746436612, + -699.6350170215012, + -76.41590901456709, + 649.7643483239231, + 145.0148701142897, + 152.66555899298714, + 65.93603488544677, + ] + ) + actual = _asm1_rhs(state, feed, kla=180.0, volume=1333.0) + np.testing.assert_allclose(actual, expected, rtol=1e-12, atol=1e-9) + + def test_settler_derivative_matches_bsm2_python_oracle(self) -> None: + components = np.array( + [30.0, 0.9, 1100.0, 50.0, 2500.0, 150.0, 450.0, + 2.0, 12.0, 1.0, 0.7, 3.5, 4.0] + ) + feed = _stream(components, 36892.0, 15.0) + soluble_indices = [0, 1, 7, 8, 9, 10, 12] + state = np.concatenate( + [ + *(feed[index] * np.linspace(0.95, 1.05, 10) + for index in soluble_indices), + np.array( + [15.0, 20.0, 30.0, 70.0, 350.0, + 500.0, 800.0, 1500.0, 3000.0, 6500.0] + ), + ] + ) + expected = np.array( + [ + 10.033888888888924, 10.033888888888924, 10.033888888888782, + 10.033888888889066, 10.247777777777287, -10.461666666666503, + -10.461666666666645, -10.461666666666503, -10.461666666666929, + -10.461666666666645, 0.3010166666666647, 0.30101666666666915, + 0.30101666666666915, 0.3010166666666647, 0.3074333333333268, + -0.3138500000000022, -0.31384999999999774, -0.31384999999999774, + -0.3138500000000066, -0.31384999999999774, 0.6689259259259295, + 0.6689259259259295, 0.6689259259259206, 0.6689259259259295, + 0.6831851851851756, -0.6974444444444394, -0.6974444444444394, + -0.6974444444444483, -0.6974444444444572, -0.6974444444444394, + 4.013555555555541, 4.013555555555612, 4.013555555555541, + 4.013555555555612, 4.099111111111, -4.184666666666672, + -4.184666666666601, -4.184666666666672, -4.184666666666743, + -4.184666666666672, 0.33446296296296474, 0.33446296296296474, + 0.3344629629629603, 0.33446296296296474, 0.3415925925925878, + -0.3487222222222197, -0.3487222222222197, -0.34872222222222415, + -0.3487222222222286, -0.3487222222222197, 0.23412407407407354, + 0.23412407407407354, 0.23412407407407354, 0.23412407407407798, + 0.23911481481481367, -0.2441055555555538, -0.2441055555555538, + -0.2441055555555538, -0.24410555555556268, -0.2441055555555538, + 1.337851851851859, 1.337851851851859, 1.3378518518518412, + 1.337851851851859, 1.3663703703703511, -1.3948888888888789, + -1.3948888888888789, -1.3948888888888966, -1.3948888888889144, + -1.3948888888888789, -159.27776512990226, -63.52879722181932, + 103.10310727314942, -476.22393432199857, 327.35724044926883, + -121215.44270969188, -208086.76214135683, -155444.6724199591, + 403397.7099670557, 73152.46245290338, + ] + ) + actual = _settler_rhs(state, feed) + np.testing.assert_allclose(actual, expected, rtol=1e-12, atol=1e-9) + def test_reference_energy_matches_published_bsm1_operating_point(self) -> None: energy = action_energy(BASELINE_ACTION) self.assertAlmostEqual(energy["aeration"], 3341.3866666666663, places=9) diff --git a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py index a3661a8e..9f6cdee7 100644 --- a/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py +++ b/benchmarks/WastewaterTreatment/BSM1AerationControl/verification/test_policy_runtime.py @@ -30,6 +30,21 @@ def test_stdout_is_isolated_from_the_json_protocol(self) -> None: runtime.reset({"scenario_id": "dry"}) self.assertEqual(runtime.control({"step": 0}), {"ok": 1}) + def test_candidate_module_supports_postponed_dataclass_annotations(self) -> None: + candidate = self._candidate( + "from __future__ import annotations\n" + "from dataclasses import dataclass\n" + "@dataclass\n" + "class State:\n" + " previous: State | None = None\n" + "state = State()\n" + "def reset_controller(scenario): state.previous = None\n" + "def control(observation): return {'loaded': state.previous is None}\n" + ) + with PolicyRuntime(candidate) as runtime: + runtime.reset({"scenario_id": "dry"}) + self.assertEqual(runtime.control({"step": 0}), {"loaded": True}) + def test_slow_control_is_terminated(self) -> None: candidate = self._candidate( "import time\n"