diff --git a/benchmarks/paper_problem1_reference.py b/benchmarks/paper_problem1_reference.py new file mode 100644 index 0000000..3c9cc56 --- /dev/null +++ b/benchmarks/paper_problem1_reference.py @@ -0,0 +1,328 @@ +# Copyright (C) 2026, SECQUOIA + +"""Generate and compare Paper Problem 1 upstream reference trajectories.""" + +from __future__ import annotations + +import argparse +import json +import os +import shlex +import subprocess +import tempfile +from pathlib import Path +from typing import Sequence + +from lyopronto.pyomo_models.paper_ocp import ( + PaperDiscretization, + classify_paper_policies, + compare_paper_problem1_trajectories, + load_upstream_matlab_trajectory, + solve_paper_problem1, +) + +DEFAULT_UPSTREAM_ROOT = Path( + os.environ.get( + "LYOPRONTO_UPSTREAM_ROOT", + "/home/bernalde/repos/simDAE-optimalcontrol-lyo", + ) +) +DEFAULT_OUTPUT = Path("benchmarks/results/paper_problem1_upstream_reference.mat") +MATLAB_RUNNER = "run_paper_problem1_upstream_reference" + + +def write_matlab_batch_files(work_dir: Path, upstream_root: Path) -> dict[str, Path]: + """Write batch-safe MATLAB helper files for the upstream OCP workflow.""" + work_dir.mkdir(parents=True, exist_ok=True) + files = { + "runner": work_dir / f"{MATLAB_RUNNER}.m", + "max_t": work_dir / "SimPy_MaxT.m", + "max_flux": work_dir / "SimPy_MaxFlux.m", + } + files["runner"].write_text(_runner_source(), encoding="utf-8") + files["max_t"].write_text( + _python_wrapper_source("SimPy_MaxT", "pyfun_MaxT.py", upstream_root), + encoding="utf-8", + ) + files["max_flux"].write_text( + _python_wrapper_source("SimPy_MaxFlux", "pyfun_MaxFlux.py", upstream_root), + encoding="utf-8", + ) + return files + + +def build_matlab_batch_command( + matlab_cmd: str, + work_dir: Path, + upstream_root: Path, + output: Path, + case_name: str = "Case2", +) -> list[str]: + """Return the MATLAB command that runs the generated batch wrapper.""" + batch = ( + f"addpath({_matlab_string(work_dir)}, '-begin'); " + f"{MATLAB_RUNNER}(" + f"{_matlab_string(upstream_root)}, " + f"{_matlab_string(output)}, " + f"{_matlab_string(case_name)});" + ) + return [matlab_cmd, "-batch", batch] + + +def compare_pyomo_to_upstream( + upstream_mat: Path, + discretization: PaperDiscretization, + solver_options: dict[str, float | int], + require_success: bool = True, +) -> dict[str, float | None]: + """Solve Paper Problem 1 from an upstream warm start and compare trajectories.""" + upstream_trajectory = load_upstream_matlab_trajectory(upstream_mat) + pyomo_result = solve_paper_problem1( + discretization=discretization, + initialization=upstream_trajectory, + solver_options=solver_options, + require_success=require_success, + ) + if "policies" not in pyomo_result: + pyomo_result["policies"] = classify_paper_policies(pyomo_result) + return compare_paper_problem1_trajectories(pyomo_result, upstream_trajectory) + + +def _runner_source() -> str: + return """function run_paper_problem1_upstream_reference(upstream_root, output_path, case_name) +code_root = fullfile(upstream_root, 'Code (Conference Version)'); +runner_root = fileparts(mfilename('fullpath')); +addpath(runner_root, '-begin'); +addpath(fullfile(code_root, 'Input Data')); +addpath(fullfile(code_root, 'Model Equations')); +addpath(fullfile(code_root, 'Events')); +addpath(fullfile(code_root, 'Simulations')); +addpath(fullfile(code_root, 'Calculations')); +addpath(fullfile(code_root, 'Sim_DAE')); + +old_dir = pwd; +cleanup = onCleanup(@() cd(old_dir)); +cd(code_root); + +ip0 = get_inputdata; +ip0 = overwrite_inputdata(ip0, string(case_name)); +ip = input_processing(ip0); +sol = Sim_1stDrying_OCP(ip); + +t = sol.t; +T = sol.T; +S = sol.S; +Tb = sol.Tb; +dSdt = sol.dSdt; +policy = sol.policy; +tsw = sol.tsw; + +output_dir = fileparts(output_path); +if ~isempty(output_dir) && ~exist(output_dir, 'dir') + mkdir(output_dir); +end +save(output_path, 't', 'T', 'S', 'Tb', 'dSdt', 'policy', 'tsw', '-v7'); +end +""" + + +def _python_wrapper_source( + function_name: str, + python_file: str, + upstream_root: Path, +) -> str: + python_dir = upstream_root / "Code (Conference Version)" / "Python" + return f"""function [Tb_opt, t_opt, S_opt] = {function_name}(ip_py) +py_file = {_matlab_string(python_file)}; +python_dir = {_matlab_string(python_dir)}; + +old_dir = pwd; +cleanup = onCleanup(@() cd(old_dir)); +cd(python_dir); + +output_py = pyrunfile(py_file, 'output_MATLAB', lyo_mpc=ip_py); +Tb_opt = double(output_py{{1}})'; +t_opt = double(output_py{{2}})'; +S_opt = double(output_py{{3}})'; + +if numel(Tb_opt) >= 3 && t_opt(3) ~= t_opt(2) + slope = (Tb_opt(3) - Tb_opt(2)) / (t_opt(3) - t_opt(2)); + Tb_opt(1) = Tb_opt(2) - slope * (t_opt(2) - t_opt(1)); +end +end +""" + + +def _matlab_string(value: str | Path) -> str: + return "'" + str(value).replace("'", "''") + "'" + + +def _validate_upstream_root(upstream_root: Path) -> None: + required = ( + upstream_root + / "Code (Conference Version)" + / "Simulations" + / "Sim_1stDrying_OCP.m" + ) + if not required.is_file(): + raise FileNotFoundError( + f"upstream root does not contain expected file: {required}" + ) + + +def _trajectory_summary(mat_path: Path) -> dict[str, float | int | None]: + trajectory = load_upstream_matlab_trajectory(mat_path) + switch_times = trajectory.get("policies", {}).get("switch_times_hr", []) + switch_times = list(switch_times) + return { + "n_points": int(len(trajectory["states"]["time_s"])), + "drying_time_hr": float(trajectory["metrics"]["drying_time_hr"]), + "first_switch_time_hr": ( + float(switch_times[0]) if len(switch_times) > 0 else None + ), + "terminal_interface_position_m": float( + trajectory["metrics"]["terminal_interface_position_m"] + ), + "max_product_temperature_K": float( + trajectory["metrics"]["max_product_temperature_K"] + ), + } + + +def _solver_options(args: argparse.Namespace) -> dict[str, float | int]: + options: dict[str, float | int] = { + "print_level": args.print_level, + "max_iter": args.max_iter, + } + if args.tol is not None: + options["tol"] = args.tol + if args.acceptable_tol is not None: + options["acceptable_tol"] = args.acceptable_tol + if args.acceptable_iter is not None: + options["acceptable_iter"] = args.acceptable_iter + return options + + +def _add_pyomo_solve_options(parser: argparse.ArgumentParser) -> None: + parser.add_argument("--n-z", type=int, default=20) + parser.add_argument("--nfe", type=int, default=12) + parser.add_argument("--ncp", type=int, default=3) + parser.add_argument("--terminal-drying-fraction", type=float, default=0.995) + parser.add_argument("--max-iter", type=int, default=3000) + parser.add_argument("--tol", type=float, default=1.0e-5) + parser.add_argument("--acceptable-tol", type=float, default=1.0e-3) + parser.add_argument("--acceptable-iter", type=int, default=5) + parser.add_argument("--print-level", type=int, default=0) + parser.add_argument( + "--allow-unsuccessful-pyomo", + action="store_true", + help="Extract and compare the Pyomo trajectory even if IPOPT is not optimal.", + ) + + +def _discretization_from_args(args: argparse.Namespace) -> PaperDiscretization: + return PaperDiscretization( + n_z=args.n_z, + nfe=args.nfe, + ncp=args.ncp, + terminal_drying_fraction=args.terminal_drying_fraction, + ) + + +def _prepare_work_dir(args: argparse.Namespace) -> tuple[Path, tempfile.TemporaryDirectory | None]: + if args.work_dir is not None: + work_dir = args.work_dir + work_dir.mkdir(parents=True, exist_ok=True) + return work_dir, None + if args.runner_only or args.keep_workdir: + return Path(tempfile.mkdtemp(prefix="lyopronto-paper-problem1-")), None + temp_dir = tempfile.TemporaryDirectory(prefix="lyopronto-paper-problem1-") + return Path(temp_dir.name), temp_dir + + +def _run_generate(args: argparse.Namespace) -> int: + upstream_root = args.upstream_root.resolve() + output = args.output.resolve() + _validate_upstream_root(upstream_root) + work_dir, cleanup = _prepare_work_dir(args) + try: + write_matlab_batch_files(work_dir, upstream_root) + command = build_matlab_batch_command( + args.matlab, + work_dir, + upstream_root, + output, + args.case, + ) + print(f"MATLAB work dir: {work_dir}") + print(f"Command: {shlex.join(command)}") + if args.runner_only: + return 0 + + subprocess.run(command, check=True) + print(json.dumps(_trajectory_summary(output), indent=2, sort_keys=True)) + if args.compare_pyomo: + comparison = compare_pyomo_to_upstream( + output, + _discretization_from_args(args), + _solver_options(args), + require_success=not args.allow_unsuccessful_pyomo, + ) + print(json.dumps(comparison, indent=2, sort_keys=True)) + finally: + if cleanup is not None: + cleanup.cleanup() + return 0 + + +def _run_compare_pyomo(args: argparse.Namespace) -> int: + comparison = compare_pyomo_to_upstream( + args.upstream_mat, + _discretization_from_args(args), + _solver_options(args), + require_success=not args.allow_unsuccessful_pyomo, + ) + print(json.dumps(comparison, indent=2, sort_keys=True)) + return 0 + + +def build_arg_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + description=__doc__, + ) + subparsers = parser.add_subparsers(dest="command", required=True) + + generate = subparsers.add_parser( + "generate", + help="Run upstream MATLAB/GEKKO and export a Paper Problem 1 .mat file.", + ) + generate.add_argument("--upstream-root", type=Path, default=DEFAULT_UPSTREAM_ROOT) + generate.add_argument("--output", type=Path, default=DEFAULT_OUTPUT) + generate.add_argument("--case", default="Case2") + generate.add_argument("--matlab", default=os.environ.get("MATLAB", "matlab")) + generate.add_argument("--work-dir", type=Path, default=None) + generate.add_argument("--runner-only", action="store_true") + generate.add_argument("--keep-workdir", action="store_true") + generate.add_argument("--compare-pyomo", action="store_true") + _add_pyomo_solve_options(generate) + generate.set_defaults(func=_run_generate) + + compare = subparsers.add_parser( + "compare-pyomo", + help="Solve Pyomo from an upstream .mat warm start and print deviations.", + ) + compare.add_argument("upstream_mat", type=Path) + _add_pyomo_solve_options(compare) + compare.set_defaults(func=_run_compare_pyomo) + + return parser + + +def main(argv: Sequence[str] | None = None) -> int: + parser = build_arg_parser() + args = parser.parse_args(argv) + return args.func(args) + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/docs/PAPER_OCP_VALIDATION.md b/docs/PAPER_OCP_VALIDATION.md index 4150aa5..79166ae 100644 --- a/docs/PAPER_OCP_VALIDATION.md +++ b/docs/PAPER_OCP_VALIDATION.md @@ -87,6 +87,45 @@ segment, Pyomo shelf temperatures were within `0.26 K` at the switch and within that offset explains the small interface-position differences in the GEKKO segment comparison. +## Regenerating The Upstream Reference + +Issue #27 is handled by `benchmarks/paper_problem1_reference.py`. The generator +keeps the upstream clone read-only: it writes temporary MATLAB wrappers for +`SimPy_MaxT` and `SimPy_MaxFlux` that use the known upstream `Python/` folder +instead of `matlab.desktop.editor.getActiveFilename`, then runs +`Sim_1stDrying_OCP` for `Case2` and saves: + +- `t` +- `T` +- `S` +- `Tb` +- `dSdt` +- `policy` +- `tsw` + +The MATLAB Python environment must be able to import GEKKO because the upstream +Policy 2 segment calls `pyfun_MaxT.py`. + +```bash +python benchmarks/paper_problem1_reference.py generate \ + --upstream-root /home/bernalde/repos/simDAE-optimalcontrol-lyo \ + --output benchmarks/results/paper_problem1_upstream_reference.mat +``` + +Use `--runner-only --work-dir /tmp/lyopronto-paper-problem1` to inspect the +generated MATLAB files and exact `matlab -batch` command without running +MATLAB. + +The exported artifact can seed the Pyomo solve and report Pyomo-vs-upstream +deviations for drying time, first switch time, terminal interface position, peak +temperature, and max-temperature profile: + +```bash +python benchmarks/paper_problem1_reference.py compare-pyomo \ + benchmarks/results/paper_problem1_upstream_reference.mat \ + --n-z 20 --nfe 12 --ncp 3 +``` + ## Mesh Diagnostics All rows use `nfe=12`, `ncp=3`, `LAGRANGE-RADAU`, the policy initializer, and a @@ -101,13 +140,12 @@ All rows use `nfe=12`, `ncp=3`, `LAGRANGE-RADAU`, the policy initializer, and a Next steps are tracked in GitHub issues: -1. #27 - Make upstream reference trajectory generation reproducible. -2. #28 - Prepare the first Paper Problem 1 validation PR. -3. #29 - Add Problem 2 with the interface-velocity constraint and expected +1. #28 - Prepare the first Paper Problem 1 validation PR. +2. #29 - Add Problem 2 with the interface-velocity constraint and expected Policy 3 -> Policy 1 -> Policy 2 sequence. -4. #30 - Compare the paper-reference transcription against LyoPRONTO's existing +3. #30 - Compare the paper-reference transcription against LyoPRONTO's existing quasi-steady Pyomo and scipy optimizers. -5. #31 - If the benchmark is credible, add a LyoPRONTO-facing experimental +4. #31 - If the benchmark is credible, add a LyoPRONTO-facing experimental policy API with the same rich result format. #26 is addressed by the bottom-node temperature constraint alignment, the diff --git a/lyopronto/pyomo_models/README.md b/lyopronto/pyomo_models/README.md index 7332451..1f5573d 100644 --- a/lyopronto/pyomo_models/README.md +++ b/lyopronto/pyomo_models/README.md @@ -92,6 +92,7 @@ drying OCP in Srisuma and Braatz, arXiv:2509.10826v1. - `generate_problem1_policy_initialization()` - Build a policy-based warm start - `initialize_paper_problem1_from_trajectory()` - Seed a model from a trajectory - `load_upstream_matlab_trajectory()` - Read upstream MATLAB output saved to `.mat` +- `compare_paper_problem1_trajectories()` - Compare Pyomo and upstream metrics - `classify_paper_policies()` - Infer active Policy 1/Policy 2 regions This module uses SI/Kelvin units from the paper/upstream code and is separate @@ -101,6 +102,24 @@ by a slow validation test using IPOPT acceptable termination (`acceptable_tol=1e-3`, `acceptable_iter=5`). Both reproduce the expected Policy 1 -> Policy 2 sequence near the paper's reported switch time. +Regenerate the full upstream Problem 1 reference with: + +```bash +python benchmarks/paper_problem1_reference.py generate \ + --upstream-root /home/bernalde/repos/simDAE-optimalcontrol-lyo \ + --output benchmarks/results/paper_problem1_upstream_reference.mat +``` + +The command writes batch-safe MATLAB wrappers outside the upstream clone and +requires MATLAB Python to import GEKKO for the upstream Policy 2 segment. Compare +the exported artifact against a Pyomo solve with: + +```bash +python benchmarks/paper_problem1_reference.py compare-pyomo \ + benchmarks/results/paper_problem1_upstream_reference.mat \ + --n-z 20 --nfe 12 --ncp 3 +``` + ## Installation ### 1. Install Pyomo diff --git a/lyopronto/pyomo_models/__init__.py b/lyopronto/pyomo_models/__init__.py index 0cf8c7b..ec43570 100644 --- a/lyopronto/pyomo_models/__init__.py +++ b/lyopronto/pyomo_models/__init__.py @@ -49,6 +49,7 @@ def _is_pyomo_available() -> bool: PaperDiscretization, PaperPrimaryDryingConfig, classify_paper_policies, + compare_paper_problem1_trajectories, create_paper_problem1_model, generate_problem1_policy_initialization, initialize_paper_problem1_from_trajectory, @@ -77,6 +78,7 @@ def _is_pyomo_available() -> bool: "generate_problem1_policy_initialization", "initialize_paper_problem1_from_trajectory", "load_upstream_matlab_trajectory", + "compare_paper_problem1_trajectories", "solve_paper_problem1", "classify_paper_policies", ] diff --git a/lyopronto/pyomo_models/paper_ocp.py b/lyopronto/pyomo_models/paper_ocp.py index f54702b..447e122 100644 --- a/lyopronto/pyomo_models/paper_ocp.py +++ b/lyopronto/pyomo_models/paper_ocp.py @@ -514,9 +514,10 @@ def load_upstream_matlab_trajectory(mat_path: str | Path) -> dict[str, Any]: """Load an upstream MATLAB trajectory saved from ``Sim_1stDrying_OCP``. The loader accepts either arrays named like the upstream output fields - (``t``, ``T``, ``S``, ``Tb``) or the first-segment arrays produced by the - local verification script (``t``, ``y``, ``Tb``), where ``y`` contains - temperature columns followed by interface position. + (``t``, ``T``, ``S``, ``Tb``, optionally ``dSdt``, ``policy``, and + ``tsw``) or the first-segment arrays produced by the local verification + script (``t``, ``y``, ``Tb``), where ``y`` contains temperature columns + followed by interface position. """ from scipy.io import loadmat @@ -561,8 +562,19 @@ def load_upstream_matlab_trajectory(mat_path: str | Path) -> dict[str, Any]: ).reshape(-1) else: interface_velocity_values = np.gradient(interface_position, time_s, edge_order=1) - - return { + if len(interface_velocity_values) != len(time_s): + raise ValueError("dSdt must share the same length as t") + + config = PaperPrimaryDryingConfig() + policies: dict[str, Any] = {} + if "policy" in data: + policies["raw_policy"] = np.atleast_1d(np.asarray(data["policy"])).reshape(-1) + if "tsw" in data: + switch_times_s = np.atleast_1d(np.asarray(data["tsw"], dtype=float)).reshape(-1) + policies["switch_times_s"] = switch_times_s + policies["switch_times_hr"] = switch_times_s / 3600.0 + + trajectory = { "states": { "time_s": time_s, "time_hr": time_s / 3600.0, @@ -574,11 +586,153 @@ def load_upstream_matlab_trajectory(mat_path: str | Path) -> dict[str, Any]: "controls": { "shelf_temperature_K": shelf_temperature, }, + "metrics": { + "drying_time_s": float(time_s[-1]), + "drying_time_hr": float(time_s[-1] / 3600.0), + "terminal_interface_position_m": float(interface_position[-1]), + "max_product_temperature_K": float(temperature.max()), + }, + "problem": { + "name": "paper_problem_1_upstream_reference", + "temperature_limit_K": config.problem1_temperature_limit, + "shelf_temperature_min_K": config.shelf_temperature_min, + "shelf_temperature_max_K": config.shelf_temperature_max, + }, "metadata": { "source": "upstream_matlab", "path": str(mat_path), }, } + if policies: + trajectory["policies"] = policies + return trajectory + + +def compare_paper_problem1_trajectories( + pyomo_result: Mapping[str, Any], + upstream_trajectory: Mapping[str, Any], +) -> dict[str, Any]: + """Compare a Pyomo Paper Problem 1 result against an upstream trajectory.""" + pyomo_drying_time_hr = _drying_time_hr(pyomo_result) + upstream_drying_time_hr = _drying_time_hr(upstream_trajectory) + pyomo_switch_time_hr = _first_switch_time_hr(pyomo_result) + upstream_switch_time_hr = _first_switch_time_hr(upstream_trajectory) + pyomo_terminal_s = _terminal_interface_position(pyomo_result) + upstream_terminal_s = _terminal_interface_position(upstream_trajectory) + pyomo_peak_temperature = _peak_temperature(pyomo_result) + upstream_peak_temperature = _peak_temperature(upstream_trajectory) + + return { + "pyomo_drying_time_hr": pyomo_drying_time_hr, + "upstream_drying_time_hr": upstream_drying_time_hr, + "drying_time_deviation_hr": pyomo_drying_time_hr - upstream_drying_time_hr, + "pyomo_first_switch_time_hr": pyomo_switch_time_hr, + "upstream_first_switch_time_hr": upstream_switch_time_hr, + "first_switch_time_deviation_hr": _nullable_difference( + pyomo_switch_time_hr, + upstream_switch_time_hr, + ), + "pyomo_terminal_interface_position_m": pyomo_terminal_s, + "upstream_terminal_interface_position_m": upstream_terminal_s, + "terminal_interface_position_deviation_m": pyomo_terminal_s - upstream_terminal_s, + "pyomo_peak_temperature_K": pyomo_peak_temperature, + "upstream_peak_temperature_K": upstream_peak_temperature, + "peak_temperature_deviation_K": pyomo_peak_temperature + - upstream_peak_temperature, + "max_temperature_profile_max_abs_deviation_K": ( + _max_temperature_profile_deviation(pyomo_result, upstream_trajectory) + ), + } + + +def _drying_time_hr(result: Mapping[str, Any]) -> float: + metrics = result.get("metrics", {}) + if "drying_time_hr" in metrics: + return float(metrics["drying_time_hr"]) + return float(np.asarray(result["states"]["time_s"], dtype=float)[-1] / 3600.0) + + +def _first_switch_time_hr(result: Mapping[str, Any]) -> float | None: + policies = result.get("policies", {}) + switch_times = policies.get("switch_times_hr") + if switch_times is not None: + switch_array = np.atleast_1d(np.asarray(switch_times, dtype=float)).reshape(-1) + if len(switch_array) > 0: + return float(switch_array[0]) + + segments = policies.get("segments", []) + if len(segments) > 1: + return float(segments[1]["start_time_hr"]) + return None + + +def _terminal_interface_position(result: Mapping[str, Any]) -> float: + metrics = result.get("metrics", {}) + if "terminal_interface_position_m" in metrics: + return float(metrics["terminal_interface_position_m"]) + interface_position = np.asarray( + result["states"]["interface_position_m"], + dtype=float, + ) + return float(interface_position[-1]) + + +def _peak_temperature(result: Mapping[str, Any]) -> float: + metrics = result.get("metrics", {}) + if "max_product_temperature_K" in metrics: + return float(metrics["max_product_temperature_K"]) + max_temperature = np.asarray(result["states"]["max_temperature_K"], dtype=float) + return float(max_temperature.max()) + + +def _nullable_difference(left: float | None, right: float | None) -> float | None: + if left is None or right is None: + return None + return float(left - right) + + +def _max_temperature_profile_deviation( + pyomo_result: Mapping[str, Any], + upstream_trajectory: Mapping[str, Any], +) -> float: + pyomo_time, pyomo_max_temperature = _unique_time_series( + pyomo_result["states"]["time_s"], + pyomo_result["states"]["max_temperature_K"], + ) + upstream_time, upstream_max_temperature = _unique_time_series( + upstream_trajectory["states"]["time_s"], + upstream_trajectory["states"]["max_temperature_K"], + ) + start = max(float(pyomo_time[0]), float(upstream_time[0])) + stop = min(float(pyomo_time[-1]), float(upstream_time[-1])) + if stop <= start: + return float("nan") + + sample_time = upstream_time[(upstream_time >= start) & (upstream_time <= stop)] + sample_time = np.unique(np.concatenate(([start], sample_time, [stop]))) + if len(sample_time) < 2: + sample_time = np.linspace(start, stop, 2) + + pyomo_interp = np.interp(sample_time, pyomo_time, pyomo_max_temperature) + upstream_interp = np.interp( + sample_time, + upstream_time, + upstream_max_temperature, + ) + return float(np.max(np.abs(pyomo_interp - upstream_interp))) + + +def _unique_time_series( + time_values: Iterable[float], + values: Iterable[float], +) -> tuple[np.ndarray, np.ndarray]: + time = np.asarray(time_values, dtype=float).reshape(-1) + series = np.asarray(values, dtype=float).reshape(-1) + order = np.argsort(time) + time = time[order] + series = series[order] + unique_time, unique_indices = np.unique(time, return_index=True) + return unique_time, series[unique_indices] def create_paper_problem1_model( diff --git a/tests/test_pyomo_models/test_init.py b/tests/test_pyomo_models/test_init.py index 1b2b32e..faa101c 100644 --- a/tests/test_pyomo_models/test_init.py +++ b/tests/test_pyomo_models/test_init.py @@ -32,6 +32,7 @@ def test_pyomo_exports_optimizers_when_pyomo_available(): "generate_problem1_policy_initialization", "initialize_paper_problem1_from_trajectory", "load_upstream_matlab_trajectory", + "compare_paper_problem1_trajectories", "solve_paper_problem1", "classify_paper_policies", } diff --git a/tests/test_pyomo_models/test_paper_ocp.py b/tests/test_pyomo_models/test_paper_ocp.py index c852abb..1a30d03 100644 --- a/tests/test_pyomo_models/test_paper_ocp.py +++ b/tests/test_pyomo_models/test_paper_ocp.py @@ -8,6 +8,7 @@ PaperDiscretization, PaperPrimaryDryingConfig, classify_paper_policies, + compare_paper_problem1_trajectories, create_paper_problem1_model, derive_primary_drying_parameters, generate_problem1_policy_initialization, @@ -196,6 +197,120 @@ def test_load_upstream_matlab_trajectory_from_segment_file(tmp_path): assert np.allclose(trajectory["controls"]["shelf_temperature_K"], 273.0) +def test_load_upstream_matlab_trajectory_preserves_full_reference_fields(tmp_path): + from scipy.io import savemat + + mat_path = tmp_path / "upstream_full.mat" + t = np.array([0.0, 3600.0, 7200.0]) + temperature = np.array( + [ + [228.0, 228.5, 229.0], + [239.0, 240.0, 241.0], + [241.0, 242.0, 243.0], + ] + ) + interface_position = np.array([0.0, 1.0e-3, 2.0e-3]) + savemat( + mat_path, + { + "t": t, + "T": temperature, + "S": interface_position, + "Tb": np.array([273.0, 270.0, 266.0]), + "dSdt": np.array([2.0e-7, 3.0e-7, 4.0e-7]), + "policy": np.array([1, 2]), + "tsw": np.array([3600.0]), + }, + ) + + trajectory = load_upstream_matlab_trajectory(mat_path) + + assert trajectory["states"]["temperature_K"].shape == (3, 3) + assert np.allclose(trajectory["states"]["interface_velocity_m_per_s"], [2.0e-7, 3.0e-7, 4.0e-7]) + assert np.allclose(trajectory["policies"]["raw_policy"], [1, 2]) + assert np.allclose(trajectory["policies"]["switch_times_hr"], [1.0]) + assert trajectory["metrics"]["drying_time_hr"] == 2.0 + assert trajectory["metrics"]["terminal_interface_position_m"] == 2.0e-3 + assert trajectory["problem"]["temperature_limit_K"] == 243.0 + + +def test_compare_paper_problem1_trajectories_reports_required_metrics(): + upstream = { + "states": { + "time_s": np.array([0.0, 3600.0, 7200.0]), + "max_temperature_K": np.array([228.0, 241.0, 243.0]), + "interface_position_m": np.array([0.0, 1.0e-3, 2.0e-3]), + }, + "metrics": { + "drying_time_hr": 2.0, + "terminal_interface_position_m": 2.0e-3, + "max_product_temperature_K": 243.0, + }, + "policies": {"switch_times_hr": np.array([1.0])}, + } + pyomo_result = { + "states": { + "time_s": np.array([0.0, 3600.0, 7560.0]), + "max_temperature_K": np.array([228.0, 241.5, 243.2]), + "interface_position_m": np.array([0.0, 1.1e-3, 2.1e-3]), + }, + "metrics": { + "drying_time_hr": 2.1, + "terminal_interface_position_m": 2.1e-3, + "max_product_temperature_K": 243.2, + }, + "policies": {"switch_times_hr": [1.1]}, + } + + comparison = compare_paper_problem1_trajectories(pyomo_result, upstream) + + assert np.isclose(comparison["drying_time_deviation_hr"], 0.1) + assert np.isclose(comparison["first_switch_time_deviation_hr"], 0.1) + assert np.isclose(comparison["terminal_interface_position_deviation_m"], 1.0e-4) + assert np.isclose(comparison["peak_temperature_deviation_K"], 0.2) + assert comparison["max_temperature_profile_max_abs_deviation_K"] > 0.0 + + +def test_paper_problem1_reference_runner_writes_batch_safe_wrappers(tmp_path): + import importlib.util + from pathlib import Path + + script_path = ( + Path(__file__).resolve().parents[2] + / "benchmarks" + / "paper_problem1_reference.py" + ) + spec = importlib.util.spec_from_file_location( + "paper_problem1_reference", + script_path, + ) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + + upstream_root = tmp_path / "upstream" + work_dir = tmp_path / "work" + + files = module.write_matlab_batch_files(work_dir, upstream_root) + command = module.build_matlab_batch_command( + "matlab", + work_dir, + upstream_root, + tmp_path / "reference.mat", + ) + + wrapper_source = files["max_t"].read_text(encoding="utf-8") + runner_source = files["runner"].read_text(encoding="utf-8") + assert "matlab.desktop.editor.getActiveFilename" not in wrapper_source + assert "pyfun_MaxT.py" in wrapper_source + assert "pyrunfile" in wrapper_source + assert "Sim_1stDrying_OCP" in runner_source + assert "policy" in runner_source + assert "tsw" in runner_source + assert command[0] == "matlab" + assert command[1] == "-batch" + assert "run_paper_problem1_upstream_reference" in command[2] + + @pytest.mark.slow @pytest.mark.skipif(not _ipopt_available(), reason="IPOPT solver not available") def test_problem1_coarse_solve_reaches_terminal_target_and_classifies_policy():