diff --git a/.gitignore b/.gitignore index a9cba54..1d70edb 100644 --- a/.gitignore +++ b/.gitignore @@ -34,6 +34,9 @@ docs/build environment.yml # for local settings .vscode/ + +# Pixi local environments +.pixi/ ################################################################################ # Jupyter Notebook # ################################################################################ diff --git a/docs/PAPER_OCP_VALIDATION.md b/docs/PAPER_OCP_VALIDATION.md new file mode 100644 index 0000000..7d048dd --- /dev/null +++ b/docs/PAPER_OCP_VALIDATION.md @@ -0,0 +1,117 @@ +# Paper OCP Validation Benchmark + +This benchmark validates whether LyoPRONTO's Pyomo direct-transcription stack +can solve published lyophilization optimal-control problems. It is intentionally +experimental and separate from the production LyoPRONTO primary-drying APIs. + +## Current Milestone + +The first implemented target is Problem 1 from Srisuma and Braatz, +arXiv:2509.10826v1: + +- primary drying only; +- moving-boundary one-dimensional frozen-region model; +- shelf temperature as the optimized control; +- objective: minimize drying time; +- constraints: product temperature at or below 243 K and shelf temperature + between 228 K and 273 K; +- expected qualitative policy sequence: maximum heat input, then product + temperature tracking. + +The implementation lives in `lyopronto.pyomo_models.paper_ocp` and uses SI +units throughout. The Pyomo benchmark itself does not depend on GEKKO or +MATLAB. GEKKO is available in the repo-local Pixi environment for upstream +policy-segment verification. + +## Validation Strategy + +The benchmark has two layers: + +- Fast tests validate parameter translation, equation helpers, model structure, + collocation discretization, and policy classification. +- Slow tests run both a coarse IPOPT solve and the upstream spatial mesh + (`n_z=20`) to verify feasibility, terminal drying, path-constraint + satisfaction, and Policy 1 -> Policy 2 detection. + +The IPOPT solve is initialized by default from a policy-based trajectory: +maximum shelf temperature until the product-temperature event, then an +algebraic Policy 2 trajectory that holds the bottom product temperature at the +limit. A MATLAB `.mat` trajectory exported from the upstream implementation can +also be loaded with `load_upstream_matlab_trajectory()` and passed to +`solve_paper_problem1(initialization=trajectory)`. + +## Upstream Solve Strategy + +The upstream implementation does not solve Problem 1 as a single generic OCP +NLP. It uses a hybrid active-policy procedure: + +1. MATLAB integrates the physical primary-drying ODE with `ode15s`. +2. Event functions detect the first violated active constraint. +3. For Policy 1, the shelf temperature is fixed at `Tbmax`. +4. For Policy 2, a GEKKO DAE simulation solves for the algebraic shelf + temperature profile that keeps the bottom product temperature at `Tmax`. +5. MATLAB re-simulates the physical ODE using that shelf profile until the next + event or drying completion. + +The GEKKO segment uses `IMODE=7`, `NODES=3`, and 30 time points for the +temperature-constrained segment. This is closer to an active-set DAE +decomposition than a one-shot free-final-time transcription. + +The current default solve uses a coarse validated spatial mesh (`n_z=5`) with a +near-complete terminal drying fraction. It checks the paper's reported switch +time near 2.4 h and a first-pass drying-time target near 6.2 h. The refined +`n_z=20` spatial mesh now also terminates cleanly with IPOPT when using the +documented acceptable NLP tolerance (`acceptable_tol=1e-3`, +`acceptable_iter=5`). The trajectory-level checks remain tight: terminal drying +within `1e-7 m`, product-temperature violation within `2e-6 K`, drying time near +6.19 h, and switch time near 2.4 h. + +As verification against the upstream reference implementation +(`PrakitrSrisuma/simDAE-optimalcontrol-lyo` commit +`5bcfece23128be7e5be51b73693dc6674223ccc6`), the MATLAB Policy 1 segment for +Problem 1 (`Case2`) was reproduced by running +`Code (Conference Version)/Simulations/Sim_1stDrying_OCP.m`. It detected the +product-temperature event at `2.363310733077 h`, which is what the policy +initializer checks against. + +The upstream GEKKO Policy 2 script (`pyfun_MaxT.py`) was also executed through +`pixi run python` with inputs built from the paper config and the policy +initializer's switch state. The local Pixi environment resolved GEKKO `1.3.2` +and successfully solved a short temperature-constrained segment with a declining +shelf-temperature profile. + +A longer GEKKO Policy 2 segment, initialized at the Problem 1 switch state and +run over the remaining Pyomo `n_z=20` drying time, gives the same shelf-control +shape as the direct-transcription result. At equal fractions of the Policy 2 +segment, Pyomo shelf temperatures were within `0.26 K` at the switch and within +`0.11 K` after the first quarter of the segment. The Pyomo switch time was +`2.3946 h` versus the policy initializer's upstream event time of `2.3633 h`; +that offset explains the small interface-position differences in the GEKKO +segment comparison. + +## Mesh Diagnostics + +All rows use `nfe=12`, `ncp=3`, `LAGRANGE-RADAU`, the policy initializer, and a +99.5% terminal drying target. + +| Spatial nodes | IPOPT status | Drying time | Policy status | +| --- | --- | --- | --- | +| `n_z=5` | optimal | ~6.19 h | Policy 1 -> Policy 2 | +| `n_z=20` | optimal/acceptable | ~6.19 h | Policy 1 -> Policy 2 | + +## Future Work + +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 + Policy 3 -> Policy 1 -> Policy 2 sequence. +4. #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 + policy API with the same rich result format. + +#26 is addressed by the bottom-node temperature constraint alignment, the +smaller expression-based NLP for vapor pressure/resistance/flux/interface +velocity, constraint scaling, and the `n_z=20` slow validation test. diff --git a/lyopronto/pyomo_models/README.md b/lyopronto/pyomo_models/README.md index 66893a2..7332451 100644 --- a/lyopronto/pyomo_models/README.md +++ b/lyopronto/pyomo_models/README.md @@ -82,6 +82,25 @@ Utility functions for initialization, scaling, and validation. - `check_solution_validity()` - Validate physical constraints - `add_scaling_suffix()` - Add scaling for numerical conditioning +### `paper_ocp.py` +Experimental paper-reference direct-transcription benchmark for the primary +drying OCP in Srisuma and Braatz, arXiv:2509.10826v1. + +**Key functions:** +- `create_paper_problem1_model()` - Build the orthogonal-collocation Pyomo model +- `solve_paper_problem1()` - Solve Paper Problem 1 with IPOPT +- `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` +- `classify_paper_policies()` - Infer active Policy 1/Policy 2 regions + +This module uses SI/Kelvin units from the paper/upstream code and is separate +from LyoPRONTO's cm/Torr/degC production APIs. The validated default solve uses +a coarse `n_z=5` mesh, and the upstream paper's `n_z=20` spatial mesh is covered +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. + ## Installation ### 1. Install Pyomo diff --git a/lyopronto/pyomo_models/__init__.py b/lyopronto/pyomo_models/__init__.py index a209f24..0cf8c7b 100644 --- a/lyopronto/pyomo_models/__init__.py +++ b/lyopronto/pyomo_models/__init__.py @@ -45,6 +45,16 @@ def _is_pyomo_available() -> bool: optimize_Pch_Tsh_pyomo, optimize_Tsh_pyomo, ) + from .paper_ocp import ( + PaperDiscretization, + PaperPrimaryDryingConfig, + classify_paper_policies, + create_paper_problem1_model, + generate_problem1_policy_initialization, + initialize_paper_problem1_from_trajectory, + load_upstream_matlab_trajectory, + solve_paper_problem1, + ) from .single_step import ( create_single_step_model, optimize_single_step, @@ -61,6 +71,14 @@ def _is_pyomo_available() -> bool: "optimize_Tsh_pyomo", "optimize_Pch_pyomo", "optimize_Pch_Tsh_pyomo", + "PaperPrimaryDryingConfig", + "PaperDiscretization", + "create_paper_problem1_model", + "generate_problem1_policy_initialization", + "initialize_paper_problem1_from_trajectory", + "load_upstream_matlab_trajectory", + "solve_paper_problem1", + "classify_paper_policies", ] __version__ = "0.1.0-dev" diff --git a/lyopronto/pyomo_models/paper_ocp.py b/lyopronto/pyomo_models/paper_ocp.py new file mode 100644 index 0000000..faea0cc --- /dev/null +++ b/lyopronto/pyomo_models/paper_ocp.py @@ -0,0 +1,1336 @@ +# Copyright (C) 2026, SECQUOIA + +"""Paper-reference optimal-control benchmarks for lyophilization. + +This module is intentionally separate from the production LyoPRONTO primary +drying models. It translates the primary-drying benchmark from Srisuma and +Braatz, arXiv:2509.10826v1, into a Pyomo direct-transcription problem so we can +validate the Pyomo.DAE + orthogonal collocation + IPOPT stack on a published +lyophilization OCP. + +The equations and defaults in this file use the paper/upstream SI convention: +temperatures in K, pressure in Pa, length in m, and time in s. Do not mix these +helpers with the existing LyoPRONTO cm/Torr/degC/cal APIs without an explicit +adapter. +""" + +from __future__ import annotations + +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Any, Iterable, Mapping + +import numpy as np + + +@dataclass(frozen=True) +class PaperPrimaryDryingConfig: + """Default primary-drying parameters from the paper's upstream code. + + The defaults are translated from: + - ``Code (Conference Version)/Input Data/get_inputdata.m`` + - ``Code (Conference Version)/Input Data/input_processing.m`` + + The first milestone targets Paper Problem 1: + minimize drying time subject to ``T(z,t) <= 243 K`` and + ``228 K <= Tb(t) <= 273 K``. + """ + + # Material properties. + solute_mass_fraction: float = 0.05 + solute_density: float = 1587.9 + water_density: float = 1000.0 + ice_density: float = 917.0 + dried_region_density: float = 215.0 + solute_heat_capacity: float = 1240.0 + ice_heat_capacity: float = 2108.0 + solute_conductivity: float = 0.126 + ice_conductivity: float = 2.25 + heat_of_sublimation: float = 2.84e6 + + # Geometry and radiation. + frozen_volume: float = 3.0e-6 + vial_diameter: float = 0.024 + glass_emissivity: float = 0.8 + steel_emissivity: float = 0.3 + stefan_boltzmann: float = 5.67e-8 + side_transfer_factor_multiplier: float = 0.78 + top_transfer_factor_multiplier: float = 1.0 + + # Primary drying conditions. + bottom_heat_transfer_coefficient: float = 25.0 + initial_temperature: float = 228.0 + wall_temperature: float = 240.0 + top_surface_temperature: float = 240.0 + initial_interface_position: float = 0.0 + chamber_water_pressure: float = 3.0 + resistance_0: float = 1.5e4 + resistance_1: float = 3.0e7 + resistance_2: float = 1.0e1 + vapor_pressure_a: float = -6139.9 + vapor_pressure_b: float = 28.8912 + microwave_heat_input: float = 0.0 + + # Paper Problem 1 OCP settings. + shelf_temperature_min: float = 228.0 + shelf_temperature_max: float = 273.0 + problem1_temperature_limit: float = 243.0 + problem1_time_guess: float = 7.0 * 3600.0 + time_bounds: tuple[float, float] = (0.25 * 3600.0, 15.0 * 3600.0) + + +@dataclass(frozen=True) +class PaperDiscretization: + """Discretization controls for the paper-reference Pyomo model.""" + + n_z: int = 5 + nfe: int = 12 + ncp: int = 3 + terminal_drying_fraction: float = 0.995 + temperature_lower_bound: float = 220.0 + temperature_upper_bound: float = 280.0 + scheme: str = "LAGRANGE-RADAU" + + +@dataclass(frozen=True) +class PaperDerivedParameters: + """Derived parameters used by the paper-reference primary drying model.""" + + solution_density: float + frozen_density: float + frozen_heat_capacity: float + frozen_conductivity: float + frozen_diffusivity: float + initial_mass: float + frozen_volume: float + cross_section_area: float + product_height: float + side_area: float + side_transfer_factor: float + top_transfer_factor: float + dpsi: float + psi: tuple[float, ...] + + +def derive_primary_drying_parameters( + config: PaperPrimaryDryingConfig, + n_z: int, +) -> PaperDerivedParameters: + """Return the upstream-derived primary-drying parameters. + + This mirrors the relevant pieces of upstream ``input_processing.m``. + """ + if n_z < 3: + raise ValueError("n_z must be at least 3 for the MOL stencil") + + xs = config.solute_mass_fraction + solution_density = 1.0 / ( + xs / config.solute_density + (1.0 - xs) / config.water_density + ) + frozen_density = 1.0 / ( + xs / config.solute_density + (1.0 - xs) / config.ice_density + ) + frozen_heat_capacity = ( + xs * config.solute_heat_capacity + + (1.0 - xs) * config.ice_heat_capacity + ) + frozen_conductivity = ( + xs * config.solute_conductivity + (1.0 - xs) * config.ice_conductivity + ) + frozen_diffusivity = frozen_conductivity / ( + frozen_density * frozen_heat_capacity + ) + initial_mass = config.frozen_volume * solution_density + frozen_volume = initial_mass / frozen_density + cross_section_area = np.pi * config.vial_diameter**2 / 4.0 + product_height = frozen_volume / cross_section_area + side_area = np.pi * config.vial_diameter * product_height + side_transfer_factor = ( + config.side_transfer_factor_multiplier * config.glass_emissivity + ) + top_transfer_factor = ( + config.top_transfer_factor_multiplier * config.glass_emissivity + ) + dpsi = 1.0 / (n_z - 1) + psi = tuple(float(i * dpsi) for i in range(n_z)) + + return PaperDerivedParameters( + solution_density=solution_density, + frozen_density=frozen_density, + frozen_heat_capacity=frozen_heat_capacity, + frozen_conductivity=frozen_conductivity, + frozen_diffusivity=frozen_diffusivity, + initial_mass=initial_mass, + frozen_volume=frozen_volume, + cross_section_area=cross_section_area, + product_height=product_height, + side_area=side_area, + side_transfer_factor=side_transfer_factor, + top_transfer_factor=top_transfer_factor, + dpsi=dpsi, + psi=psi, + ) + + +def saturation_pressure( + temperature: np.ndarray | float, + config: PaperPrimaryDryingConfig = PaperPrimaryDryingConfig(), +) -> np.ndarray | float: + """Return water saturation pressure [Pa] from temperature [K].""" + return np.exp(config.vapor_pressure_a / temperature + config.vapor_pressure_b) + + +def product_resistance( + interface_position: np.ndarray | float, + config: PaperPrimaryDryingConfig = PaperPrimaryDryingConfig(), +) -> np.ndarray | float: + """Return cake resistance [m/s] at interface position ``S`` [m].""" + return config.resistance_0 + config.resistance_1 * interface_position / ( + 1.0 + config.resistance_2 * interface_position + ) + + +def sublimation_flux( + interface_temperature: np.ndarray | float, + interface_position: np.ndarray | float, + config: PaperPrimaryDryingConfig = PaperPrimaryDryingConfig(), +) -> np.ndarray | float: + """Return sublimation flux [kg/(m^2 s)]. + + This helper mirrors the smooth signed equation used in the + direct-transcription model. The Pyomo model adds a separate constraint that + keeps the sublimation driving force nonnegative. + """ + pressure = saturation_pressure(interface_temperature, config) + resistance = product_resistance(interface_position, config) + return (pressure - config.chamber_water_pressure) / resistance + + +def interface_velocity( + interface_temperature: np.ndarray | float, + interface_position: np.ndarray | float, + config: PaperPrimaryDryingConfig = PaperPrimaryDryingConfig(), + derived: PaperDerivedParameters | None = None, +) -> np.ndarray | float: + """Return interface velocity ``dS/dt`` [m/s].""" + if derived is None: + derived = derive_primary_drying_parameters(config, n_z=20) + flux = sublimation_flux(interface_temperature, interface_position, config) + return flux / (derived.frozen_density - config.dried_region_density) + + +def generate_problem1_policy_initialization( + config: PaperPrimaryDryingConfig | None = None, + discretization: PaperDiscretization | None = None, + n_time_points: int = 240, + rtol: float = 1.0e-7, + atol: float = 1.0e-9, +) -> dict[str, Any]: + """Generate an upstream-style policy trajectory for Problem 1 warm starts. + + This mirrors the paper implementation's control logic for Problem 1 without + depending on MATLAB or GEKKO: integrate Policy 1 with maximum shelf + temperature until the bottom product temperature reaches the limit, then use + the Policy 2 algebraic condition to keep the bottom temperature at the limit. + """ + from scipy.integrate import solve_ivp + + config = config or PaperPrimaryDryingConfig() + discretization = discretization or PaperDiscretization() + derived = derive_primary_drying_parameters(config, discretization.n_z) + n_z = discretization.n_z + target_s = discretization.terminal_drying_fraction * derived.product_height + y0 = np.concatenate( + ( + np.full(n_z, config.initial_temperature, dtype=float), + np.array([config.initial_interface_position], dtype=float), + ) + ) + + def policy1_rhs(_time, y): + temperature = y[:n_z] + interface_position = float(y[-1]) + dtemperature_dt, dinterface_dt = _numeric_temperature_rhs( + temperature, + interface_position, + config.shelf_temperature_max, + config, + derived, + ) + return np.concatenate((dtemperature_dt, [dinterface_dt])) + + def temperature_event(_time, y): + return y[n_z - 1] - config.problem1_temperature_limit + + temperature_event.terminal = True + temperature_event.direction = 1 + + def finish_event(_time, y): + return y[-1] - target_s + + finish_event.terminal = True + finish_event.direction = 1 + + sol1 = solve_ivp( + policy1_rhs, + (0.0, config.time_bounds[1]), + y0, + method="BDF", + dense_output=True, + events=(temperature_event, finish_event), + rtol=rtol, + atol=atol, + max_step=100.0, + ) + if not sol1.success: + raise RuntimeError(f"Policy 1 initialization failed: {sol1.message}") + + switch_time = ( + float(sol1.t_events[0][0]) if len(sol1.t_events[0]) else float(sol1.t[-1]) + ) + finished_in_policy1 = bool(len(sol1.t_events[1])) + + sol2 = None + final_time = float(sol1.t_events[1][0]) if finished_in_policy1 else None + if not finished_in_policy1: + switch_state = sol1.sol(switch_time) + y2_0 = np.concatenate( + ( + switch_state[: n_z - 1], + np.array([float(switch_state[-1])], dtype=float), + ) + ) + + def policy2_rhs(_time, y): + temperature = np.empty(n_z, dtype=float) + temperature[: n_z - 1] = y[: n_z - 1] + temperature[n_z - 1] = config.problem1_temperature_limit + interface_position = float(y[-1]) + shelf_temperature = _policy2_shelf_temperature( + temperature, + interface_position, + config, + derived, + ) + shelf_temperature = float( + np.clip( + shelf_temperature, + config.shelf_temperature_min, + config.shelf_temperature_max, + ) + ) + dtemperature_dt, dinterface_dt = _numeric_temperature_rhs( + temperature, + interface_position, + shelf_temperature, + config, + derived, + ) + return np.concatenate((dtemperature_dt[: n_z - 1], [dinterface_dt])) + + def finish_event_policy2(_time, y): + return y[-1] - target_s + + finish_event_policy2.terminal = True + finish_event_policy2.direction = 1 + + sol2 = solve_ivp( + policy2_rhs, + (switch_time, config.time_bounds[1]), + y2_0, + method="BDF", + dense_output=True, + events=finish_event_policy2, + rtol=rtol, + atol=atol, + max_step=100.0, + ) + if not sol2.success: + raise RuntimeError(f"Policy 2 initialization failed: {sol2.message}") + final_time = ( + float(sol2.t_events[0][0]) if len(sol2.t_events[0]) else float(sol2.t[-1]) + ) + + time_s = _policy_sample_times(switch_time, final_time, n_time_points, sol2 is not None) + temperature_values = np.zeros((len(time_s), n_z), dtype=float) + interface_values = np.zeros(len(time_s), dtype=float) + shelf_values = np.zeros(len(time_s), dtype=float) + interface_velocity_values = np.zeros(len(time_s), dtype=float) + flux_values = np.zeros(len(time_s), dtype=float) + labels: list[str] = [] + + for index, time in enumerate(time_s): + if sol2 is None or time <= switch_time: + state = sol1.sol(time) + temperature = np.asarray(state[:n_z], dtype=float) + interface_position = float(state[-1]) + shelf_temperature = config.shelf_temperature_max + label = "policy_1_max_heat_input" + else: + state = sol2.sol(time) + temperature = np.empty(n_z, dtype=float) + temperature[: n_z - 1] = np.asarray(state[: n_z - 1], dtype=float) + temperature[n_z - 1] = config.problem1_temperature_limit + interface_position = float(state[-1]) + shelf_temperature = _policy2_shelf_temperature( + temperature, + interface_position, + config, + derived, + ) + shelf_temperature = float( + np.clip( + shelf_temperature, + config.shelf_temperature_min, + config.shelf_temperature_max, + ) + ) + label = "policy_2_temperature_tracking" + + _, dinterface_dt = _numeric_temperature_rhs( + temperature, + interface_position, + shelf_temperature, + config, + derived, + ) + flux = float(sublimation_flux(temperature[0], interface_position, config)) + + temperature_values[index, :] = temperature + interface_values[index] = interface_position + shelf_values[index] = shelf_temperature + interface_velocity_values[index] = dinterface_dt + flux_values[index] = flux + labels.append(label) + + return { + "states": { + "time_s": time_s, + "time_hr": time_s / 3600.0, + "temperature_K": temperature_values, + "max_temperature_K": temperature_values.max(axis=1), + "interface_position_m": interface_values, + "interface_velocity_m_per_s": interface_velocity_values, + "sublimation_flux_kg_m2_s": flux_values, + }, + "controls": { + "shelf_temperature_K": shelf_values, + }, + "policies": { + "labels": labels, + "segments": _compress_policy_labels(time_s / 3600.0, labels), + "switch_times_hr": [switch_time / 3600.0] if sol2 is not None else [], + }, + "metrics": { + "drying_time_s": final_time, + "drying_time_hr": final_time / 3600.0, + "policy1_switch_time_s": switch_time, + "policy1_switch_time_hr": switch_time / 3600.0, + "terminal_drying_fraction": float(interface_values[-1] / derived.product_height), + }, + "metadata": { + "source": "problem1_policy_initialization", + "n_z": n_z, + "rtol": rtol, + "atol": atol, + }, + "problem": { + "name": "paper_problem_1_policy_initialization", + "temperature_limit_K": config.problem1_temperature_limit, + "shelf_temperature_min_K": config.shelf_temperature_min, + "shelf_temperature_max_K": config.shelf_temperature_max, + "terminal_drying_fraction_target": discretization.terminal_drying_fraction, + }, + } + + +def initialize_paper_problem1_from_trajectory( + model: Any, + trajectory: Mapping[str, Any], + set_final_time: bool = True, +) -> None: + """Initialize a discretized Paper Problem 1 model from a trajectory dict.""" + states = trajectory.get("states", trajectory) + controls = trajectory.get("controls", {}) + time_s = np.asarray(states["time_s"], dtype=float) + temperature = np.asarray(states["temperature_K"], dtype=float) + interface_position = np.asarray(states["interface_position_m"], dtype=float) + if "shelf_temperature_K" in controls: + shelf_temperature = np.asarray(controls["shelf_temperature_K"], dtype=float) + else: + shelf_temperature = np.asarray(states["shelf_temperature_K"], dtype=float) + + if temperature.ndim != 2: + raise ValueError("temperature_K must be a 2-D array with shape (time, z)") + if len(time_s) != len(temperature): + raise ValueError("time_s and temperature_K have inconsistent lengths") + if not np.all(np.diff(time_s) > 0.0): + raise ValueError("time_s must be strictly increasing") + + final_time = float(time_s[-1]) + if set_final_time: + _set_var_value_within_bounds(model.t_final, final_time) + + source_psi = np.linspace(0.0, 1.0, temperature.shape[1]) + target_psi = np.linspace(0.0, 1.0, len(list(model.z))) + dtemperature_dt = np.gradient(temperature, time_s, axis=0, edge_order=1) + dinterface_dt = ( + np.asarray(states["interface_velocity_m_per_s"], dtype=float) + if "interface_velocity_m_per_s" in states + else np.gradient(interface_position, time_s, edge_order=1) + ) + + for tau in sorted(model.t): + absolute_time = float(tau) * final_time + interface_value = float(np.interp(absolute_time, time_s, interface_position)) + shelf_value = float(np.interp(absolute_time, time_s, shelf_temperature)) + source_temperature_at_time = np.array( + [np.interp(absolute_time, time_s, temperature[:, j]) for j in range(temperature.shape[1])] + ) + source_dtemperature_at_time = np.array( + [np.interp(absolute_time, time_s, dtemperature_dt[:, j]) for j in range(temperature.shape[1])] + ) + target_temperature = np.interp(target_psi, source_psi, source_temperature_at_time) + target_dtemperature_dt = np.interp( + target_psi, + source_psi, + source_dtemperature_at_time, + ) + + _set_var_value_within_bounds(model.S[tau], interface_value) + _set_var_value_within_bounds(model.Tb[tau], shelf_value) + for i, value in zip(model.z, target_temperature): + _set_var_value_within_bounds(model.T[i, tau], float(value)) + if hasattr(model, "dT_dtau"): + model.dT_dtau[i, tau].set_value(float(final_time * target_dtemperature_dt[i])) + + velocity = float(np.interp(absolute_time, time_s, dinterface_dt)) + velocity = max(velocity, 1.0e-12) + + if hasattr(model, "dS_dtau"): + model.dS_dtau[tau].set_value(final_time * velocity) + + +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. + """ + from scipy.io import loadmat + + data = loadmat(mat_path, squeeze_me=True) + if "t" not in data: + raise ValueError("MATLAB trajectory must contain a 't' time array") + + time_s = np.atleast_1d(np.asarray(data["t"], dtype=float)).reshape(-1) + if "T" in data: + temperature = np.asarray(data["T"], dtype=float) + elif "y" in data: + y = np.asarray(data["y"], dtype=float) + if y.ndim != 2 or y.shape[1] < 2: + raise ValueError("'y' must have temperature columns followed by S") + temperature = y[:, :-1] + else: + raise ValueError("MATLAB trajectory must contain either 'T' or 'y'") + + if temperature.ndim == 1: + temperature = temperature.reshape(-1, 1) + if "S" in data: + interface_position = np.atleast_1d(np.asarray(data["S"], dtype=float)).reshape(-1) + elif "y" in data: + interface_position = np.asarray(data["y"], dtype=float)[:, -1] + else: + raise ValueError("MATLAB trajectory must contain 'S' when 'y' is absent") + + if "Tb" not in data: + raise ValueError("MATLAB trajectory must contain a 'Tb' shelf-temperature array") + shelf_temperature = np.atleast_1d(np.asarray(data["Tb"], dtype=float)).reshape(-1) + + if len(time_s) != len(temperature) or len(time_s) != len(interface_position): + raise ValueError("MATLAB trajectory arrays must share the same time dimension") + if len(shelf_temperature) == 1: + shelf_temperature = np.full_like(time_s, float(shelf_temperature[0])) + if len(shelf_temperature) != len(time_s): + raise ValueError("Tb must be scalar or share the same length as t") + + if "dSdt" in data: + interface_velocity_values = np.atleast_1d( + np.asarray(data["dSdt"], dtype=float) + ).reshape(-1) + else: + interface_velocity_values = np.gradient(interface_position, time_s, edge_order=1) + + return { + "states": { + "time_s": time_s, + "time_hr": time_s / 3600.0, + "temperature_K": temperature, + "max_temperature_K": temperature.max(axis=1), + "interface_position_m": interface_position, + "interface_velocity_m_per_s": interface_velocity_values, + }, + "controls": { + "shelf_temperature_K": shelf_temperature, + }, + "metadata": { + "source": "upstream_matlab", + "path": str(mat_path), + }, + } + + +def create_paper_problem1_model( + config: PaperPrimaryDryingConfig | None = None, + discretization: PaperDiscretization | None = None, + apply_scaling: bool = False, +): + """Create the Pyomo direct-transcription model for Paper Problem 1. + + The model is discretized immediately with orthogonal collocation so callers + receive an NLP ready for initialization and solve. Scaling suffixes are + optional because the validated coarse benchmark solves cleanly without them. + """ + import pyomo.dae as dae + import pyomo.environ as pyo + + config = config or PaperPrimaryDryingConfig() + discretization = discretization or PaperDiscretization() + derived = derive_primary_drying_parameters(config, discretization.n_z) + if not 0.0 < discretization.terminal_drying_fraction < 1.0: + raise ValueError("terminal_drying_fraction must be in (0, 1)") + + model = pyo.ConcreteModel() + model._paper_config = config + model._paper_discretization = discretization + model._paper_derived = derived + + model.z = pyo.RangeSet(0, discretization.n_z - 1) + model.t = dae.ContinuousSet(bounds=(0.0, 1.0)) + model.psi = pyo.Param( + model.z, + initialize={i: derived.psi[i] for i in range(discretization.n_z)}, + ) + + terminal_s = discretization.terminal_drying_fraction * derived.product_height + model.t_final = pyo.Var( + bounds=config.time_bounds, + initialize=config.problem1_time_guess, + ) + model.T = pyo.Var( + model.z, + model.t, + bounds=( + discretization.temperature_lower_bound, + discretization.temperature_upper_bound, + ), + initialize=config.initial_temperature, + ) + model.S = pyo.Var( + model.t, + bounds=(0.0, terminal_s), + initialize=0.5 * terminal_s, + ) + model.Tb = pyo.Var( + model.t, + bounds=(config.shelf_temperature_min, config.shelf_temperature_max), + initialize=config.shelf_temperature_max, + ) + + model.dT_dtau = dae.DerivativeVar(model.T, wrt=model.t) + model.dS_dtau = dae.DerivativeVar(model.S, wrt=model.t) + + def vapor_pressure_rule(m, t): + return pyo.exp( + config.vapor_pressure_a / m.T[0, t] + config.vapor_pressure_b + ) + + model.Pw = pyo.Expression(model.t, rule=vapor_pressure_rule) + + def resistance_rule(m, t): + return config.resistance_0 + config.resistance_1 * m.S[t] / ( + 1.0 + config.resistance_2 * m.S[t] + ) + + model.Rp = pyo.Expression(model.t, rule=resistance_rule) + + def sublimation_flux_rule(m, t): + return (m.Pw[t] - config.chamber_water_pressure) / m.Rp[t] + + model.Nw = pyo.Expression(model.t, rule=sublimation_flux_rule) + + def nonnegative_sublimation_flux_rule(m, t): + return m.Pw[t] >= config.chamber_water_pressure + + model.nonnegative_sublimation_flux = pyo.Constraint( + model.t, + rule=nonnegative_sublimation_flux_rule, + ) + + def interface_velocity_rule(m, t): + return m.Nw[t] / (derived.frozen_density - config.dried_region_density) + + model.dSdt = pyo.Expression(model.t, rule=interface_velocity_rule) + + def interface_ode_rule(m, t): + if t == m.t.first(): + return pyo.Constraint.Skip + return m.dS_dtau[t] == m.t_final * m.dSdt[t] + + model.interface_ode = pyo.Constraint(model.t, rule=interface_ode_rule) + + def temperature_rhs(m, i, t): + thickness = derived.product_height - m.S[t] + volume = derived.cross_section_area * thickness + side_loss = ( + derived.side_transfer_factor + * config.stefan_boltzmann + * derived.side_area + * (m.T[i, t] ** 4 - config.wall_temperature**4) + / (volume * derived.frozen_density * derived.frozen_heat_capacity) + ) + source = config.microwave_heat_input / ( + volume * derived.frozen_density * derived.frozen_heat_capacity + ) + + if i == 0: + top_radiation = ( + derived.top_transfer_factor + * config.stefan_boltzmann + * (m.T[i, t] ** 4 - config.top_surface_temperature**4) + ) + diffusion = ( + derived.frozen_diffusivity + / thickness**2 + / derived.dpsi**2 + * ( + 2.0 * m.T[1, t] + - 2.0 * m.T[0, t] + - 2.0 + * m.Nw[t] + * derived.dpsi + * config.heat_of_sublimation + * thickness + / derived.frozen_conductivity + - top_radiation + * 2.0 + * derived.dpsi + * thickness + / derived.frozen_conductivity + ) + ) + convection_gradient = ( + thickness + * m.Nw[t] + * config.heat_of_sublimation + / derived.frozen_conductivity + + top_radiation + * thickness + / derived.frozen_conductivity + ) + convection = ( + -((m.psi[i] - 1.0) * m.dSdt[t] / thickness) * convection_gradient + ) + return diffusion + convection - side_loss + source + + if i == discretization.n_z - 1: + diffusion = ( + derived.frozen_diffusivity + / thickness**2 + / derived.dpsi**2 + * ( + 2.0 * m.T[i - 1, t] + - 2.0 * m.T[i, t] + + 2.0 + * (m.S[t] - derived.product_height) + * config.bottom_heat_transfer_coefficient + * derived.dpsi + * (m.T[i, t] - m.Tb[t]) + / derived.frozen_conductivity + ) + ) + convection_gradient = ( + (m.S[t] - derived.product_height) + * config.bottom_heat_transfer_coefficient + * (m.T[i, t] - m.Tb[t]) + / derived.frozen_conductivity + ) + convection = ( + -((m.psi[i] - 1.0) * m.dSdt[t] / thickness) * convection_gradient + ) + return diffusion + convection - side_loss + source + + diffusion = ( + derived.frozen_diffusivity + / thickness**2 + / derived.dpsi**2 + * (m.T[i - 1, t] - 2.0 * m.T[i, t] + m.T[i + 1, t]) + ) + convection = ( + -((m.psi[i] - 1.0) * m.dSdt[t] / thickness) + * (m.T[i + 1, t] - m.T[i - 1, t]) + / (2.0 * derived.dpsi) + ) + return diffusion + convection - side_loss + source + + def temperature_ode_rule(m, i, t): + if t == m.t.first(): + return pyo.Constraint.Skip + return m.dT_dtau[i, t] == m.t_final * temperature_rhs(m, i, t) + + model.temperature_ode = pyo.Constraint( + model.z, + model.t, + rule=temperature_ode_rule, + ) + + model.initial_interface = pyo.Constraint( + expr=model.S[0.0] == config.initial_interface_position + ) + + def initial_temperature_rule(m, i): + return m.T[i, 0.0] == config.initial_temperature + + model.initial_temperature = pyo.Constraint( + model.z, + rule=initial_temperature_rule, + ) + + bottom_node = discretization.n_z - 1 + + def product_temperature_limit_rule(m, t): + return m.T[bottom_node, t] <= config.problem1_temperature_limit + + model.product_temperature_limit = pyo.Constraint( + model.t, + rule=product_temperature_limit_rule, + ) + + model.terminal_drying = pyo.Constraint(expr=model.S[1.0] >= terminal_s) + model.objective = pyo.Objective(expr=model.t_final, sense=pyo.minimize) + + discretizer = pyo.TransformationFactory("dae.collocation") + discretizer.apply_to( + model, + nfe=discretization.nfe, + ncp=discretization.ncp, + scheme=discretization.scheme, + ) + + _initialize_problem1_model(model) + if apply_scaling: + _add_problem1_scaling(model) + + return model + + +def _initialize_problem1_model(model: Any) -> None: + """Populate deterministic initial guesses after collocation.""" + config = model._paper_config + discretization = model._paper_discretization + derived = model._paper_derived + terminal_s = discretization.terminal_drying_fraction * derived.product_height + time_guess = config.problem1_time_guess + model.t_final.set_value(time_guess) + + for t in sorted(model.t): + tau = float(t) + s_guess = terminal_s * tau + model.S[t].set_value(s_guess) + + if tau <= 0.35: + tb_guess = config.shelf_temperature_max + else: + decline = (tau - 0.35) / 0.65 + tb_guess = config.shelf_temperature_max - 20.0 * min(max(decline, 0.0), 1.0) + model.Tb[t].set_value(max(config.shelf_temperature_min, tb_guess)) + + top_guess = min( + config.problem1_temperature_limit - 1.0, + config.initial_temperature + 13.0 * tau, + ) + bottom_guess = min( + config.problem1_temperature_limit - 0.25, + config.initial_temperature + 15.0 * tau, + ) + for i in model.z: + frac = i / (discretization.n_z - 1) + model.T[i, t].set_value(top_guess + frac * (bottom_guess - top_guess)) + + top_temperature = model.T[0, t].value + pressure = float(saturation_pressure(top_temperature, config)) + resistance = float(product_resistance(s_guess, config)) + flux = max((pressure - config.chamber_water_pressure) / resistance, 1.0e-8) + velocity = flux / (derived.frozen_density - config.dried_region_density) + if hasattr(model, "dS_dtau"): + model.dS_dtau[t].set_value(time_guess * velocity) + + +def _add_problem1_scaling(model: Any) -> None: + """Add IPOPT scaling suffixes for the SI-unit benchmark model.""" + import pyomo.environ as pyo + + model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + model.scaling_factor[model.t_final] = 1.0e-4 + for t in model.t: + model.scaling_factor[model.S[t]] = 1.0e2 + model.scaling_factor[model.Tb[t]] = 1.0e-2 + for i in model.z: + model.scaling_factor[model.T[i, t]] = 1.0e-2 + + _set_component_scaling(model, "interface_ode", 1.0e2) + _set_component_scaling(model, "temperature_ode", 1.0e-2) + _set_component_scaling(model, "initial_interface", 1.0e2) + _set_component_scaling(model, "initial_temperature", 1.0e-2) + _set_component_scaling(model, "product_temperature_limit", 1.0e-2) + _set_component_scaling(model, "nonnegative_sublimation_flux", 1.0e-1) + _set_component_scaling(model, "terminal_drying", 1.0e2) + _set_component_scaling(model, "dS_dtau_disc_eq", 1.0e2) + _set_component_scaling(model, "dT_dtau_disc_eq", 1.0e-2) + + +def _set_component_scaling(model: Any, component_name: str, factor: float) -> None: + if not hasattr(model, component_name): + return + component = getattr(model, component_name) + if component.is_indexed(): + for index in component: + model.scaling_factor[component[index]] = factor + else: + model.scaling_factor[component] = factor + + +def _numeric_temperature_rhs( + temperature: np.ndarray, + interface_position: float, + shelf_temperature: float, + config: PaperPrimaryDryingConfig, + derived: PaperDerivedParameters, +) -> tuple[np.ndarray, float]: + """Return dimensional RHS values for the paper primary-drying model.""" + n_z = len(temperature) + thickness = max(derived.product_height - interface_position, 1.0e-9) + volume = derived.cross_section_area * thickness + pressure = float(saturation_pressure(temperature[0], config)) + resistance = float(product_resistance(interface_position, config)) + flux = max((pressure - config.chamber_water_pressure) / resistance, 0.0) + dinterface_dt = flux / (derived.frozen_density - config.dried_region_density) + dtemperature_dt = np.zeros(n_z, dtype=float) + + for i in range(n_z): + side_loss = ( + derived.side_transfer_factor + * config.stefan_boltzmann + * derived.side_area + * (temperature[i] ** 4 - config.wall_temperature**4) + / (volume * derived.frozen_density * derived.frozen_heat_capacity) + ) + source = config.microwave_heat_input / ( + volume * derived.frozen_density * derived.frozen_heat_capacity + ) + + if i == 0: + top_radiation = ( + derived.top_transfer_factor + * config.stefan_boltzmann + * (temperature[i] ** 4 - config.top_surface_temperature**4) + ) + diffusion = ( + derived.frozen_diffusivity + / thickness**2 + / derived.dpsi**2 + * ( + 2.0 * temperature[1] + - 2.0 * temperature[0] + - 2.0 + * flux + * derived.dpsi + * config.heat_of_sublimation + * thickness + / derived.frozen_conductivity + - top_radiation + * 2.0 + * derived.dpsi + * thickness + / derived.frozen_conductivity + ) + ) + convection_gradient = ( + thickness * flux * config.heat_of_sublimation / derived.frozen_conductivity + + top_radiation * thickness / derived.frozen_conductivity + ) + convection = ( + -((derived.psi[i] - 1.0) * dinterface_dt / thickness) + * convection_gradient + ) + dtemperature_dt[i] = diffusion + convection - side_loss + source + elif i == n_z - 1: + diffusion = ( + derived.frozen_diffusivity + / thickness**2 + / derived.dpsi**2 + * ( + 2.0 * temperature[i - 1] + - 2.0 * temperature[i] + + 2.0 + * (interface_position - derived.product_height) + * config.bottom_heat_transfer_coefficient + * derived.dpsi + * (temperature[i] - shelf_temperature) + / derived.frozen_conductivity + ) + ) + convection_gradient = ( + (interface_position - derived.product_height) + * config.bottom_heat_transfer_coefficient + * (temperature[i] - shelf_temperature) + / derived.frozen_conductivity + ) + convection = ( + -((derived.psi[i] - 1.0) * dinterface_dt / thickness) + * convection_gradient + ) + dtemperature_dt[i] = diffusion + convection - side_loss + source + else: + diffusion = ( + derived.frozen_diffusivity + / thickness**2 + / derived.dpsi**2 + * ( + temperature[i - 1] + - 2.0 * temperature[i] + + temperature[i + 1] + ) + ) + convection = ( + -((derived.psi[i] - 1.0) * dinterface_dt / thickness) + * (temperature[i + 1] - temperature[i - 1]) + / (2.0 * derived.dpsi) + ) + dtemperature_dt[i] = diffusion + convection - side_loss + source + + return dtemperature_dt, dinterface_dt + + +def _policy2_shelf_temperature( + temperature: np.ndarray, + interface_position: float, + config: PaperPrimaryDryingConfig, + derived: PaperDerivedParameters, +) -> float: + """Return shelf temperature that makes the bottom-node derivative zero.""" + bottom_index = len(temperature) - 1 + thickness = max(derived.product_height - interface_position, 1.0e-9) + volume = derived.cross_section_area * thickness + bottom_temperature = temperature[bottom_index] + side_loss = ( + derived.side_transfer_factor + * config.stefan_boltzmann + * derived.side_area + * (bottom_temperature**4 - config.wall_temperature**4) + / (volume * derived.frozen_density * derived.frozen_heat_capacity) + ) + source = config.microwave_heat_input / ( + volume * derived.frozen_density * derived.frozen_heat_capacity + ) + diffusion_scale = derived.frozen_diffusivity / thickness**2 / derived.dpsi**2 + target_stencil = (side_loss - source) / diffusion_scale + bottom_delta = ( + target_stencil + - 2.0 * temperature[bottom_index - 1] + + 2.0 * bottom_temperature + ) + denominator = ( + 2.0 + * (interface_position - derived.product_height) + * config.bottom_heat_transfer_coefficient + * derived.dpsi + ) + bottom_minus_shelf = bottom_delta * derived.frozen_conductivity / denominator + return bottom_temperature - bottom_minus_shelf + + +def _policy_sample_times( + switch_time: float, + final_time: float, + n_time_points: int, + has_policy2: bool, +) -> np.ndarray: + """Return monotonically increasing sample times preserving the switch point.""" + if n_time_points < 4: + raise ValueError("n_time_points must be at least 4") + if not has_policy2: + return np.linspace(0.0, final_time, n_time_points) + + policy1_fraction = max(0.1, min(0.9, switch_time / final_time)) + n_policy1 = max(3, int(round(n_time_points * policy1_fraction))) + n_policy2 = max(3, n_time_points - n_policy1 + 1) + return np.unique( + np.concatenate( + ( + np.linspace(0.0, switch_time, n_policy1), + np.linspace(switch_time, final_time, n_policy2), + ) + ) + ) + + +def _set_var_value_within_bounds(var: Any, value: float) -> None: + """Set a Pyomo VarData value after clipping to its finite bounds.""" + lower, upper = var.bounds + if lower is not None: + value = max(float(lower), value) + if upper is not None: + value = min(float(upper), value) + var.set_value(value) + + +def solve_paper_problem1( + config: PaperPrimaryDryingConfig | None = None, + discretization: PaperDiscretization | None = None, + solver: str = "ipopt", + solver_options: Mapping[str, Any] | None = None, + initialization: str | Mapping[str, Any] | None = "policy", + tee: bool = False, + require_success: bool = True, + return_model: bool = False, +) -> dict[str, Any]: + """Build and solve Paper Problem 1 with Pyomo/IPOPT.""" + import pyomo.environ as pyo + + model = create_paper_problem1_model(config, discretization) + config = config or model._paper_config + discretization = discretization or model._paper_discretization + if initialization == "policy": + trajectory = generate_problem1_policy_initialization(config, discretization) + initialize_paper_problem1_from_trajectory(model, trajectory) + elif isinstance(initialization, Mapping): + initialize_paper_problem1_from_trajectory(model, initialization) + elif initialization is not None: + raise ValueError("initialization must be 'policy', a trajectory mapping, or None") + + try: + from idaes.core.solvers import get_solver + + opt = get_solver(solver) + except Exception: + opt = pyo.SolverFactory(solver) + + if solver == "ipopt" and hasattr(opt, "options"): + opt.options.setdefault("max_iter", 5000) + opt.options.setdefault("tol", 1.0e-6) + opt.options.setdefault("acceptable_tol", 1.0e-3) + opt.options.setdefault("acceptable_iter", 5) + opt.options.setdefault("constr_viol_tol", 1.0e-6) + opt.options.setdefault("mu_strategy", "adaptive") + opt.options.setdefault("bound_relax_factor", 1.0e-8) + opt.options.setdefault("nlp_scaling_method", "user-scaling") + opt.options.setdefault("print_level", 5 if tee else 0) + + if solver_options: + for key, value in solver_options.items(): + opt.options[key] = value + + results = opt.solve(model, tee=tee) + solution = extract_paper_solution(model, results) + solution["policies"] = classify_paper_policies(solution) + if return_model: + solution["model"] = model + + if require_success and not _is_successful_termination(results): + metadata = solution["metadata"] + raise RuntimeError( + "Paper Problem 1 solve did not converge " + f"(status={metadata['status']}, " + f"termination_condition={metadata['termination_condition']})" + ) + + return solution + + +def extract_paper_solution(model: Any, results: Any | None = None) -> dict[str, Any]: + """Extract a solved or initialized paper OCP model into a rich result dict.""" + import pyomo.environ as pyo + + config = model._paper_config + discretization = model._paper_discretization + derived = model._paper_derived + t_points = sorted(model.t) + z_points = list(model.z) + t_final = float(pyo.value(model.t_final)) + tau = np.array([float(t) for t in t_points]) + time_s = tau * t_final + + temperature = np.array( + [[float(pyo.value(model.T[i, t])) for i in z_points] for t in t_points] + ) + interface_position = np.array([float(pyo.value(model.S[t])) for t in t_points]) + shelf_temperature = np.array([float(pyo.value(model.Tb[t])) for t in t_points]) + interface_velocity_values = np.array( + [float(pyo.value(model.dSdt[t])) for t in t_points] + ) + flux = np.array([float(pyo.value(model.Nw[t])) for t in t_points]) + resistance = np.array([float(pyo.value(model.Rp[t])) for t in t_points]) + vapor_pressure = np.array([float(pyo.value(model.Pw[t])) for t in t_points]) + max_temperature = temperature.max(axis=1) + + target_s = discretization.terminal_drying_fraction * derived.product_height + metrics = { + "drying_time_s": t_final, + "drying_time_hr": t_final / 3600.0, + "terminal_interface_position_m": float(interface_position[-1]), + "terminal_drying_fraction": float(interface_position[-1] / derived.product_height), + "target_interface_position_m": target_s, + "terminal_gap_m": max(0.0, target_s - float(interface_position[-1])), + "max_product_temperature_K": float(max_temperature.max()), + "max_temperature_violation_K": max( + 0.0, + float(max_temperature.max() - config.problem1_temperature_limit), + ), + "shelf_lower_violation_K": max( + 0.0, + float(config.shelf_temperature_min - shelf_temperature.min()), + ), + "shelf_upper_violation_K": max( + 0.0, + float(shelf_temperature.max() - config.shelf_temperature_max), + ), + } + + status = None + termination_condition = None + if results is not None: + solver_info = getattr(results, "solver", None) + status = str(getattr(solver_info, "status", None)) + termination_condition = str( + getattr(solver_info, "termination_condition", None) + ) + + return { + "states": { + "tau": tau, + "time_s": time_s, + "time_hr": time_s / 3600.0, + "temperature_K": temperature, + "max_temperature_K": max_temperature, + "interface_position_m": interface_position, + "interface_velocity_m_per_s": interface_velocity_values, + "sublimation_flux_kg_m2_s": flux, + "resistance_m_per_s": resistance, + "vapor_pressure_Pa": vapor_pressure, + }, + "controls": { + "shelf_temperature_K": shelf_temperature, + }, + "metrics": metrics, + "metadata": { + "status": status, + "termination_condition": termination_condition, + "n_z": discretization.n_z, + "nfe": discretization.nfe, + "ncp": discretization.ncp, + "scheme": discretization.scheme, + }, + "problem": { + "name": "paper_problem_1", + "temperature_limit_K": config.problem1_temperature_limit, + "shelf_temperature_min_K": config.shelf_temperature_min, + "shelf_temperature_max_K": config.shelf_temperature_max, + "terminal_drying_fraction_target": discretization.terminal_drying_fraction, + }, + "config": asdict(config), + "derived": asdict(derived), + } + + +def classify_paper_policies( + result: Mapping[str, Any], + tolerances: Mapping[str, float] | None = None, +) -> dict[str, Any]: + """Infer active paper policies from a result trajectory. + + Policy 1 is maximum heat input (``Tb = Tb_max``). Policy 2 is product + temperature tracking (``max(T) = T_limit``). If both are active at a point, + Policy 2 takes precedence because it is the path constraint that forces the + control away from unconstrained maximum heating. + """ + tolerances = dict(tolerances or {}) + temperature_tolerance = tolerances.get("temperature_K", 0.35) + shelf_tolerance = tolerances.get("shelf_temperature_K", 0.35) + + states = result["states"] + controls = result["controls"] + problem = result["problem"] + time_hr = np.asarray(states["time_hr"]) + max_temperature = np.asarray(states["max_temperature_K"]) + shelf_temperature = np.asarray(controls["shelf_temperature_K"]) + temperature_limit = float(problem["temperature_limit_K"]) + shelf_max = float(problem["shelf_temperature_max_K"]) + + labels: list[str] = [] + for temp, shelf in zip(max_temperature, shelf_temperature): + temp_active = temp >= temperature_limit - temperature_tolerance + shelf_active = abs(shelf - shelf_max) <= shelf_tolerance + if temp_active: + labels.append("policy_2_temperature_tracking") + elif shelf_active: + labels.append("policy_1_max_heat_input") + else: + labels.append("unclassified") + + segments = _compress_policy_labels(time_hr, labels) + switch_times = [segment["start_time_hr"] for segment in segments[1:]] + return { + "labels": labels, + "segments": segments, + "switch_times_hr": switch_times, + } + + +def _compress_policy_labels( + time_hr: Iterable[float], + labels: Iterable[str], +) -> list[dict[str, Any]]: + """Return contiguous policy-label segments.""" + times = list(time_hr) + label_list = list(labels) + if not label_list: + return [] + + segments: list[dict[str, Any]] = [] + current = label_list[0] + start_index = 0 + for index, label in enumerate(label_list[1:], start=1): + if label != current: + segments.append( + { + "label": current, + "start_time_hr": float(times[start_index]), + "end_time_hr": float(times[index - 1]), + } + ) + current = label + start_index = index + segments.append( + { + "label": current, + "start_time_hr": float(times[start_index]), + "end_time_hr": float(times[-1]), + } + ) + return segments + + +def _is_successful_termination(results: Any) -> bool: + """Return whether a Pyomo solve status is acceptable for extraction.""" + import pyomo.environ as pyo + + solver = getattr(results, "solver", None) + termination_condition = getattr(solver, "termination_condition", None) + acceptable = {str(pyo.TerminationCondition.optimal).lower()} + locally_optimal = getattr(pyo.TerminationCondition, "locallyOptimal", None) + if locally_optimal is not None: + acceptable.add(str(locally_optimal).lower()) + return str(termination_condition).lower() in acceptable diff --git a/pixi.lock b/pixi.lock new file mode 100644 index 0000000..91a1a45 --- /dev/null +++ b/pixi.lock @@ -0,0 +1,552 @@ +version: 6 +environments: + default: + channels: + - url: https://conda.anaconda.org/conda-forge/ + indexes: + - https://pypi.org/simple + options: + pypi-prerelease-mode: if-necessary-or-explicit + packages: + linux-64: + - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2026.4.22-hbd8a1cb_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.45.1-default_hbd61a6d_102.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.8.0-hecca717_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libffi-3.5.2-h3435931_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libgcc-15.2.0-he0feb66_19.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libgomp-15.2.0-he0feb66_19.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.3-hb03c661_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libmpdec-4.0.0-hb03c661_1.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.53.1-h0c1763c_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.42-h5347b49_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.2-h25fd6f3_2.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.6-hdb14827_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/openssl-3.6.2-h35e630c_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/python-3.13.13-h6add32d_100_cp313.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.13-8_cp313.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/tk-8.6.13-noxft_h366c992_103.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2025c-hc9c84f9_1.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.7-hb78ec9c_6.conda + - pypi: https://files.pythonhosted.org/packages/4b/32/e0f13a1c5b0f8572d0ec6ae2f6c677b7991fafd95da523159c19eff0696a/contourpy-1.3.3-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/e7/05/c19819d5e3d95294a6f5947fb9b9629efb316b96de511b418c53d245aae6/cycler-0.12.1-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/e2/98/8b1e801939839d405f1f122e7d175cebe9aeb4e114f95bfc45e3152af9a7/fonttools-4.62.1-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/47/c6/4a03e540b0cbd7ea29437b547370a74732b35d23e22d20735718642a2140/gekko-1.3.2-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/2b/0a/7b98e1e119878a27ba8618ca1e18b14f992ff1eda40f47bccccf4de44121/kiwisolver-1.5.0-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/8a/17/4402d0d14ccf1dfc70932600b68097fbbf9c898a4871d2cbbe79c7801a32/matplotlib-3.10.9-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/d1/73/a9d864e42a01896bb5974475438f16086be9ba1f0d19d0bb7a07427c4a8b/numpy-2.4.4-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/df/b2/87e62e8c3e2f4b32e5fe99e0b86d576da1312593b39f47d8ceef365e95ed/packaging-26.2-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/67/ee/21d4e8536afd1a328f01b359b4d3997b291ffd35a237c877b331c1c3b71c/pillow-12.2.0-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/10/bd/c038d7cc38edc1aa5bf91ab8068b63d4308c66c4c8bb3cbba7dfbc049f9c/pyparsing-3.3.2-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/ec/57/56b9bcc3c9c6a792fcbaf139543cee77261f3651ca9da0c93f5c1221264b/python_dateutil-2.9.0.post0-py2.py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/b8/0c/51f6841f1d84f404f92463fc2b1ba0da357ca1e3db6b7fbda26956c3b82a/ruamel_yaml-0.19.1-py3-none-any.whl + - pypi: https://files.pythonhosted.org/packages/f5/5f/f17563f28ff03c7b6799c50d01d5d856a1d55f2676f537ca8d28c7f627cd/scipy-1.17.1-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + - pypi: https://files.pythonhosted.org/packages/b7/ce/149a00dd41f10bc29e5921b496af8b574d8413afcd5e30dfa0ed46c2cc5e/six-1.17.0-py2.py3-none-any.whl + - pypi: ./ +packages: +- conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda + build_number: 20 + sha256: 1dd3fffd892081df9726d7eb7e0dea6198962ba775bd88842135a4ddb4deb3c9 + md5: a9f577daf3de00bca7c3c76c0ecbd1de + depends: + - __glibc >=2.17,<3.0.a0 + - libgomp >=7.5.0 + constrains: + - openmp_impl <0.0a0 + license: BSD-3-Clause + license_family: BSD + purls: [] + size: 28948 + timestamp: 1770939786096 +- conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda + sha256: 0b75d45f0bba3e95dc693336fa51f40ea28c980131fec438afb7ce6118ed05f6 + md5: d2ffd7602c02f2b316fd921d39876885 + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + license: bzip2-1.0.6 + license_family: BSD + purls: [] + size: 260182 + timestamp: 1771350215188 +- conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2026.4.22-hbd8a1cb_0.conda + sha256: c9dbcc8039a52023660d6d1bbf87594a93dd69c6ac5a2a44323af2c92976728d + md5: e18ad67cf881dcadee8b8d9e2f8e5f73 + depends: + - __unix + license: ISC + purls: [] + size: 131039 + timestamp: 1776865545798 +- pypi: https://files.pythonhosted.org/packages/4b/32/e0f13a1c5b0f8572d0ec6ae2f6c677b7991fafd95da523159c19eff0696a/contourpy-1.3.3-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + name: contourpy + version: 1.3.3 + sha256: 4debd64f124ca62069f313a9cb86656ff087786016d76927ae2cf37846b006c9 + requires_dist: + - numpy>=1.25 + - furo ; extra == 'docs' + - sphinx>=7.2 ; extra == 'docs' + - sphinx-copybutton ; extra == 'docs' + - bokeh ; extra == 'bokeh' + - selenium ; extra == 'bokeh' + - contourpy[bokeh,docs] ; extra == 'mypy' + - bokeh ; extra == 'mypy' + - docutils-stubs ; extra == 'mypy' + - mypy==1.17.0 ; extra == 'mypy' + - types-pillow ; extra == 'mypy' + - contourpy[test-no-images] ; extra == 'test' + - matplotlib ; extra == 'test' + - pillow ; extra == 'test' + - pytest ; extra == 'test-no-images' + - pytest-cov ; extra == 'test-no-images' + - pytest-rerunfailures ; extra == 'test-no-images' + - pytest-xdist ; extra == 'test-no-images' + - wurlitzer ; extra == 'test-no-images' + requires_python: '>=3.11' +- pypi: https://files.pythonhosted.org/packages/e7/05/c19819d5e3d95294a6f5947fb9b9629efb316b96de511b418c53d245aae6/cycler-0.12.1-py3-none-any.whl + name: cycler + version: 0.12.1 + sha256: 85cef7cff222d8644161529808465972e51340599459b8ac3ccbac5a854e0d30 + requires_dist: + - ipython ; extra == 'docs' + - matplotlib ; extra == 'docs' + - numpydoc ; extra == 'docs' + - sphinx ; extra == 'docs' + - pytest ; extra == 'tests' + - pytest-cov ; extra == 'tests' + - pytest-xdist ; extra == 'tests' + requires_python: '>=3.8' +- pypi: https://files.pythonhosted.org/packages/e2/98/8b1e801939839d405f1f122e7d175cebe9aeb4e114f95bfc45e3152af9a7/fonttools-4.62.1-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl + name: fonttools + version: 4.62.1 + sha256: 6706d1cb1d5e6251a97ad3c1b9347505c5615c112e66047abbef0f8545fa30d1 + requires_dist: + - lxml>=4.0 ; extra == 'lxml' + - brotli>=1.0.1 ; platform_python_implementation == 'CPython' and extra == 'woff' + - brotlicffi>=0.8.0 ; platform_python_implementation != 'CPython' and extra == 'woff' + - zopfli>=0.1.4 ; extra == 'woff' + - unicodedata2>=17.0.0 ; python_full_version < '3.15' and extra == 'unicode' + - lz4>=1.7.4.2 ; extra == 'graphite' + - scipy ; platform_python_implementation != 'PyPy' and extra == 'interpolatable' + - munkres ; platform_python_implementation == 'PyPy' and extra == 'interpolatable' + - pycairo ; extra == 'interpolatable' + - matplotlib ; extra == 'plot' + - sympy ; extra == 'symfont' + - xattr ; sys_platform == 'darwin' and extra == 'type1' + - skia-pathops>=0.5.0 ; extra == 'pathops' + - uharfbuzz>=0.45.0 ; extra == 'repacker' + - lxml>=4.0 ; extra == 'all' + - brotli>=1.0.1 ; platform_python_implementation == 'CPython' and extra == 'all' + - brotlicffi>=0.8.0 ; platform_python_implementation != 'CPython' and extra == 'all' + - zopfli>=0.1.4 ; extra == 'all' + - unicodedata2>=17.0.0 ; python_full_version < '3.15' and extra == 'all' + - lz4>=1.7.4.2 ; extra == 'all' + - scipy ; platform_python_implementation != 'PyPy' and extra == 'all' + - munkres ; platform_python_implementation == 'PyPy' and extra == 'all' + - pycairo ; extra == 'all' + - matplotlib ; extra == 'all' + - sympy ; extra == 'all' + - xattr ; sys_platform == 'darwin' and extra == 'all' + - skia-pathops>=0.5.0 ; extra == 'all' + - uharfbuzz>=0.45.0 ; extra == 'all' + requires_python: '>=3.10' +- pypi: https://files.pythonhosted.org/packages/47/c6/4a03e540b0cbd7ea29437b547370a74732b35d23e22d20735718642a2140/gekko-1.3.2-py3-none-any.whl + name: gekko + version: 1.3.2 + sha256: d1ae9a096eb496f16df8f42a3dba7dcea59e2d17bf9a7d6d2b9187a06959e00a + requires_dist: + - numpy>=1.8 + requires_python: '>=2.6' +- pypi: https://files.pythonhosted.org/packages/2b/0a/7b98e1e119878a27ba8618ca1e18b14f992ff1eda40f47bccccf4de44121/kiwisolver-1.5.0-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl + name: kiwisolver + version: 1.5.0 + sha256: 332b4f0145c30b5f5ad9374881133e5aa64320428a57c2c2b61e9d891a51c2f3 + requires_python: '>=3.10' +- conda: https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.45.1-default_hbd61a6d_102.conda + sha256: 3d584956604909ff5df353767f3a2a2f60e07d070b328d109f30ac40cd62df6c + md5: 18335a698559cdbcd86150a48bf54ba6 + depends: + - __glibc >=2.17,<3.0.a0 + - zstd >=1.5.7,<1.6.0a0 + constrains: + - binutils_impl_linux-64 2.45.1 + license: GPL-3.0-only + license_family: GPL + purls: [] + size: 728002 + timestamp: 1774197446916 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.8.0-hecca717_0.conda + sha256: ea33c40977ea7a2c3658c522230058395bc2ee0d89d99f0711390b6a1ee80d12 + md5: a3b390520c563d78cc58974de95a03e5 + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + constrains: + - expat 2.8.0.* + license: MIT + license_family: MIT + purls: [] + size: 77241 + timestamp: 1777846112704 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libffi-3.5.2-h3435931_0.conda + sha256: 31f19b6a88ce40ebc0d5a992c131f57d919f73c0b92cd1617a5bec83f6e961e6 + md5: a360c33a5abe61c07959e449fa1453eb + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + license: MIT + license_family: MIT + purls: [] + size: 58592 + timestamp: 1769456073053 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libgcc-15.2.0-he0feb66_19.conda + sha256: 8e0a3b5e41272e5678499b5dfc4cddb673f9e935de01eb0767ce857001229f46 + md5: 57736f29cc2b0ec0b6c2952d3f101b6a + depends: + - __glibc >=2.17,<3.0.a0 + - _openmp_mutex >=4.5 + constrains: + - libgcc-ng ==15.2.0=*_19 + - libgomp 15.2.0 he0feb66_19 + license: GPL-3.0-only WITH GCC-exception-3.1 + license_family: GPL + purls: [] + size: 1041084 + timestamp: 1778269013026 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libgomp-15.2.0-he0feb66_19.conda + sha256: 5abe4ab9d93f6c9757d654f1969ae2267d4505315c1f2f8fe705fd60af084f1b + md5: faac990cb7aedc7f3a2224f2c9b0c26c + depends: + - __glibc >=2.17,<3.0.a0 + license: GPL-3.0-only WITH GCC-exception-3.1 + license_family: GPL + purls: [] + size: 603817 + timestamp: 1778268942614 +- conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.3-hb03c661_0.conda + sha256: ec30e52a3c1bf7d0425380a189d209a52baa03f22fb66dd3eb587acaa765bd6d + md5: b88d90cad08e6bc8ad540cb310a761fb + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + constrains: + - xz 5.8.3.* + license: 0BSD + purls: [] + size: 113478 + timestamp: 1775825492909 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libmpdec-4.0.0-hb03c661_1.conda + sha256: fe171ed5cf5959993d43ff72de7596e8ac2853e9021dec0344e583734f1e0843 + md5: 2c21e66f50753a083cbe6b80f38268fa + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + license: BSD-2-Clause + license_family: BSD + purls: [] + size: 92400 + timestamp: 1769482286018 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.53.1-h0c1763c_0.conda + sha256: 54cdcd3214313b62c2a8ee277e6f42150d9b748264c1b70d958bf735e420ef8d + md5: 7dc38adcbf71e6b38748e919e16e0dce + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + - libzlib >=1.3.2,<2.0a0 + license: blessing + purls: [] + size: 954962 + timestamp: 1777986471789 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.42-h5347b49_0.conda + sha256: bc1b08c92626c91500fd9f26f2c797f3eb153b627d53e9c13cd167f1e12b2829 + md5: 38ffe67b78c9d4de527be8315e5ada2c + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + license: BSD-3-Clause + license_family: BSD + purls: [] + size: 40297 + timestamp: 1775052476770 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.2-h25fd6f3_2.conda + sha256: 55044c403570f0dc26e6364de4dc5368e5f3fc7ff103e867c487e2b5ab2bcda9 + md5: d87ff7921124eccd67248aa483c23fec + depends: + - __glibc >=2.17,<3.0.a0 + constrains: + - zlib 1.3.2 *_2 + license: Zlib + license_family: Other + purls: [] + size: 63629 + timestamp: 1774072609062 +- pypi: ./ + name: lyopronto + version: 1.1.0 + sha256: e65f5c547d4ff8072a852109c003796b4c9a6bf23a67436ba69a4b48155b719f + requires_dist: + - numpy>=1.24.0 + - scipy>=1.10.0 + - matplotlib>=3.7.0 + - ruamel-yaml>=0.18.0 + - pyomo>=6.7.0 ; extra == 'optimization' + - idaes-pse>=2.5,<2.6 ; python_full_version < '3.10' and extra == 'optimization' + - idaes-pse>=2.9.0 ; python_full_version >= '3.10' and extra == 'optimization' + - pytest>=7.4.0 ; extra == 'dev' + - pytest-mock>=3 ; extra == 'dev' + - pytest-cov>=4.1.0 ; extra == 'dev' + - pytest-xdist>=3.3.0 ; extra == 'dev' + - hypothesis>=6.82.0 ; extra == 'dev' + - ruff>=0.12.0 ; extra == 'dev' + - mypy>=1.4.0 ; extra == 'dev' + - pandas>=2.0 ; extra == 'dev' + - papermill>=2.6.0 ; extra == 'dev' + - ipykernel>=6.15.0 ; extra == 'dev' + - mkdocstrings-python>=2 ; extra == 'docs' + - mkdocs-material>=9.1.0 ; extra == 'docs' + - mike>=2.1 ; extra == 'docs' + - mkdocs-ipynb>=0.1.1 ; extra == 'docs' + - lyopronto[optimization,dev,docs] ; extra == 'all' + requires_python: '>=3.8' +- pypi: https://files.pythonhosted.org/packages/8a/17/4402d0d14ccf1dfc70932600b68097fbbf9c898a4871d2cbbe79c7801a32/matplotlib-3.10.9-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl + name: matplotlib + version: 3.10.9 + sha256: 8f3bcac1ca5ed000a6f4337d47ba67dfddf37ed6a46c15fd7f014997f7bf865f + requires_dist: + - contourpy>=1.0.1 + - cycler>=0.10 + - fonttools>=4.22.0 + - kiwisolver>=1.3.1 + - numpy>=1.23 + - packaging>=20.0 + - pillow>=8 + - pyparsing>=3 + - python-dateutil>=2.7 + - meson-python>=0.13.1,<0.17.0 ; extra == 'dev' + - pybind11>=2.13.2,!=2.13.3 ; extra == 'dev' + - setuptools-scm>=7,<10 ; extra == 'dev' + - setuptools>=64 ; extra == 'dev' + requires_python: '>=3.10' +- conda: https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.6-hdb14827_0.conda + sha256: fc89f74bbe362fb29fa3c037697a89bec140b346a2469a90f7936d1d7ea4d8a3 + md5: fc21868a1a5aacc937e7a18747acb8a5 + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + license: X11 AND BSD-3-Clause + purls: [] + size: 918956 + timestamp: 1777422145199 +- pypi: https://files.pythonhosted.org/packages/d1/73/a9d864e42a01896bb5974475438f16086be9ba1f0d19d0bb7a07427c4a8b/numpy-2.4.4-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + name: numpy + version: 2.4.4 + sha256: c901b15172510173f5cb310eae652908340f8dede90fff9e3bf6c0d8dfd92f83 + requires_python: '>=3.11' +- conda: https://conda.anaconda.org/conda-forge/linux-64/openssl-3.6.2-h35e630c_0.conda + sha256: c0ef482280e38c71a08ad6d71448194b719630345b0c9c60744a2010e8a8e0cb + md5: da1b85b6a87e141f5140bb9924cecab0 + depends: + - __glibc >=2.17,<3.0.a0 + - ca-certificates + - libgcc >=14 + license: Apache-2.0 + license_family: Apache + purls: [] + size: 3167099 + timestamp: 1775587756857 +- pypi: https://files.pythonhosted.org/packages/df/b2/87e62e8c3e2f4b32e5fe99e0b86d576da1312593b39f47d8ceef365e95ed/packaging-26.2-py3-none-any.whl + name: packaging + version: '26.2' + sha256: 5fc45236b9446107ff2415ce77c807cee2862cb6fac22b8a73826d0693b0980e + requires_python: '>=3.8' +- pypi: https://files.pythonhosted.org/packages/67/ee/21d4e8536afd1a328f01b359b4d3997b291ffd35a237c877b331c1c3b71c/pillow-12.2.0-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + name: pillow + version: 12.2.0 + sha256: eedf4b74eda2b5a4b2b2fb4c006d6295df3bf29e459e198c90ea48e130dc75c3 + requires_dist: + - furo ; extra == 'docs' + - olefile ; extra == 'docs' + - sphinx>=8.2 ; extra == 'docs' + - sphinx-autobuild ; extra == 'docs' + - sphinx-copybutton ; extra == 'docs' + - sphinx-inline-tabs ; extra == 'docs' + - sphinxext-opengraph ; extra == 'docs' + - olefile ; extra == 'fpx' + - olefile ; extra == 'mic' + - arro3-compute ; extra == 'test-arrow' + - arro3-core ; extra == 'test-arrow' + - nanoarrow ; extra == 'test-arrow' + - pyarrow ; extra == 'test-arrow' + - check-manifest ; extra == 'tests' + - coverage>=7.4.2 ; extra == 'tests' + - defusedxml ; extra == 'tests' + - markdown2 ; extra == 'tests' + - olefile ; extra == 'tests' + - packaging ; extra == 'tests' + - pyroma>=5 ; extra == 'tests' + - pytest ; extra == 'tests' + - pytest-cov ; extra == 'tests' + - pytest-timeout ; extra == 'tests' + - pytest-xdist ; extra == 'tests' + - trove-classifiers>=2024.10.12 ; extra == 'tests' + - defusedxml ; extra == 'xmp' + requires_python: '>=3.10' +- pypi: https://files.pythonhosted.org/packages/10/bd/c038d7cc38edc1aa5bf91ab8068b63d4308c66c4c8bb3cbba7dfbc049f9c/pyparsing-3.3.2-py3-none-any.whl + name: pyparsing + version: 3.3.2 + sha256: 850ba148bd908d7e2411587e247a1e4f0327839c40e2e5e6d05a007ecc69911d + requires_dist: + - railroad-diagrams ; extra == 'diagrams' + - jinja2 ; extra == 'diagrams' + requires_python: '>=3.9' +- conda: https://conda.anaconda.org/conda-forge/linux-64/python-3.13.13-h6add32d_100_cp313.conda + build_number: 100 + sha256: 7f77eb57648f545c1f58e10035d0d9d66b0a0efb7c4b58d3ed89ec7269afdde1 + md5: 05051be49267378d2fcd12931e319ac3 + depends: + - __glibc >=2.17,<3.0.a0 + - bzip2 >=1.0.8,<2.0a0 + - ld_impl_linux-64 >=2.36.1 + - libexpat >=2.7.5,<3.0a0 + - libffi >=3.5.2,<3.6.0a0 + - libgcc >=14 + - liblzma >=5.8.2,<6.0a0 + - libmpdec >=4.0.0,<5.0a0 + - libsqlite >=3.52.0,<4.0a0 + - libuuid >=2.42,<3.0a0 + - libzlib >=1.3.2,<2.0a0 + - ncurses >=6.5,<7.0a0 + - openssl >=3.5.6,<4.0a0 + - python_abi 3.13.* *_cp313 + - readline >=8.3,<9.0a0 + - tk >=8.6.13,<8.7.0a0 + - tzdata + license: Python-2.0 + purls: [] + size: 37358322 + timestamp: 1775614712638 + python_site_packages_path: lib/python3.13/site-packages +- pypi: https://files.pythonhosted.org/packages/ec/57/56b9bcc3c9c6a792fcbaf139543cee77261f3651ca9da0c93f5c1221264b/python_dateutil-2.9.0.post0-py2.py3-none-any.whl + name: python-dateutil + version: 2.9.0.post0 + sha256: a8b2bc7bffae282281c8140a97d3aa9c14da0b136dfe83f850eea9a5f7470427 + requires_dist: + - six>=1.5 + requires_python: '>=2.7,!=3.0.*,!=3.1.*,!=3.2.*' +- conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.13-8_cp313.conda + build_number: 8 + sha256: 210bffe7b121e651419cb196a2a63687b087497595c9be9d20ebe97dd06060a7 + md5: 94305520c52a4aa3f6c2b1ff6008d9f8 + constrains: + - python 3.13.* *_cp313 + license: BSD-3-Clause + license_family: BSD + purls: [] + size: 7002 + timestamp: 1752805902938 +- conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda + sha256: 12ffde5a6f958e285aa22c191ca01bbd3d6e710aa852e00618fa6ddc59149002 + md5: d7d95fc8287ea7bf33e0e7116d2b95ec + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + - ncurses >=6.5,<7.0a0 + license: GPL-3.0-only + license_family: GPL + purls: [] + size: 345073 + timestamp: 1765813471974 +- pypi: https://files.pythonhosted.org/packages/b8/0c/51f6841f1d84f404f92463fc2b1ba0da357ca1e3db6b7fbda26956c3b82a/ruamel_yaml-0.19.1-py3-none-any.whl + name: ruamel-yaml + version: 0.19.1 + sha256: 27592957fedf6e0b62f281e96effd28043345e0e66001f97683aa9a40c667c93 + requires_dist: + - ruamel-yaml-clib ; platform_python_implementation == 'CPython' and extra == 'oldlibyaml' + - ruamel-yaml-clibz>=0.3.7 ; platform_python_implementation == 'CPython' and extra == 'libyaml' + - ruamel-yaml-jinja2>=0.2 ; extra == 'jinja2' + - ryd ; extra == 'docs' + - mercurial>5.7 ; extra == 'docs' + requires_python: '>=3.9' +- pypi: https://files.pythonhosted.org/packages/f5/5f/f17563f28ff03c7b6799c50d01d5d856a1d55f2676f537ca8d28c7f627cd/scipy-1.17.1-cp313-cp313-manylinux_2_27_x86_64.manylinux_2_28_x86_64.whl + name: scipy + version: 1.17.1 + sha256: 581b2264fc0aa555f3f435a5944da7504ea3a065d7029ad60e7c3d1ae09c5464 + requires_dist: + - numpy>=1.26.4,<2.7 + - pytest>=8.0.0 ; extra == 'test' + - pytest-cov ; extra == 'test' + - pytest-timeout ; extra == 'test' + - pytest-xdist ; extra == 'test' + - asv ; extra == 'test' + - mpmath ; extra == 'test' + - gmpy2 ; extra == 'test' + - threadpoolctl ; extra == 'test' + - scikit-umfpack ; extra == 'test' + - pooch ; extra == 'test' + - hypothesis>=6.30 ; extra == 'test' + - array-api-strict>=2.3.1 ; extra == 'test' + - cython ; extra == 'test' + - meson ; extra == 'test' + - ninja ; sys_platform != 'emscripten' and extra == 'test' + - sphinx>=5.0.0,<8.2.0 ; extra == 'doc' + - intersphinx-registry ; extra == 'doc' + - pydata-sphinx-theme>=0.15.2 ; extra == 'doc' + - sphinx-copybutton ; extra == 'doc' + - sphinx-design>=0.4.0 ; extra == 'doc' + - matplotlib>=3.5 ; extra == 'doc' + - numpydoc ; extra == 'doc' + - jupytext ; extra == 'doc' + - myst-nb>=1.2.0 ; extra == 'doc' + - pooch ; extra == 'doc' + - jupyterlite-sphinx>=0.19.1 ; extra == 'doc' + - jupyterlite-pyodide-kernel ; extra == 'doc' + - linkify-it-py ; extra == 'doc' + - tabulate ; extra == 'doc' + - click<8.3.0 ; extra == 'dev' + - spin ; extra == 'dev' + - mypy==1.10.0 ; extra == 'dev' + - typing-extensions ; extra == 'dev' + - types-psutil ; extra == 'dev' + - pycodestyle ; extra == 'dev' + - ruff>=0.12.0 ; extra == 'dev' + - cython-lint>=0.12.2 ; extra == 'dev' + requires_python: '>=3.11' +- pypi: https://files.pythonhosted.org/packages/b7/ce/149a00dd41f10bc29e5921b496af8b574d8413afcd5e30dfa0ed46c2cc5e/six-1.17.0-py2.py3-none-any.whl + name: six + version: 1.17.0 + sha256: 4721f391ed90541fddacab5acf947aa0d3dc7d27b2e1e8eda2be8970586c3274 + requires_python: '>=2.7,!=3.0.*,!=3.1.*,!=3.2.*' +- conda: https://conda.anaconda.org/conda-forge/linux-64/tk-8.6.13-noxft_h366c992_103.conda + sha256: cafeec44494f842ffeca27e9c8b0c27ed714f93ac77ddadc6aaf726b5554ebac + md5: cffd3bdd58090148f4cfcd831f4b26ab + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + - libzlib >=1.3.1,<2.0a0 + constrains: + - xorg-libx11 >=1.8.12,<2.0a0 + license: TCL + license_family: BSD + purls: [] + size: 3301196 + timestamp: 1769460227866 +- conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2025c-hc9c84f9_1.conda + sha256: 1d30098909076af33a35017eed6f2953af1c769e273a0626a04722ac4acaba3c + md5: ad659d0a2b3e47e38d829aa8cad2d610 + license: LicenseRef-Public-Domain + purls: [] + size: 119135 + timestamp: 1767016325805 +- conda: https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.7-hb78ec9c_6.conda + sha256: 68f0206ca6e98fea941e5717cec780ed2873ffabc0e1ed34428c061e2c6268c7 + md5: 4a13eeac0b5c8e5b8ab496e6c4ddd829 + depends: + - __glibc >=2.17,<3.0.a0 + - libzlib >=1.3.1,<2.0a0 + license: BSD-3-Clause + license_family: BSD + purls: [] + size: 601375 + timestamp: 1764777111296 diff --git a/pyproject.toml b/pyproject.toml index 8808155..023c10f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -144,3 +144,14 @@ quote-style = "double" indent-style = "space" # Unix line endings line-ending = "lf" + +[tool.pixi.workspace] +channels = ["conda-forge"] +platforms = ["linux-64"] + +[tool.pixi.dependencies] +python = "3.13.*" + +[tool.pixi.pypi-dependencies] +lyopronto = { path = ".", editable = true } +gekko = ">=1.3.0,<2" diff --git a/tests/test_pyomo_models/test_init.py b/tests/test_pyomo_models/test_init.py index c604aa6..1b2b32e 100644 --- a/tests/test_pyomo_models/test_init.py +++ b/tests/test_pyomo_models/test_init.py @@ -26,6 +26,14 @@ def test_pyomo_exports_optimizers_when_pyomo_available(): "optimize_Tsh_pyomo", "optimize_Pch_pyomo", "optimize_Pch_Tsh_pyomo", + "PaperPrimaryDryingConfig", + "PaperDiscretization", + "create_paper_problem1_model", + "generate_problem1_policy_initialization", + "initialize_paper_problem1_from_trajectory", + "load_upstream_matlab_trajectory", + "solve_paper_problem1", + "classify_paper_policies", } if pyomo_models.PYOMO_AVAILABLE: diff --git a/tests/test_pyomo_models/test_paper_ocp.py b/tests/test_pyomo_models/test_paper_ocp.py new file mode 100644 index 0000000..a697d44 --- /dev/null +++ b/tests/test_pyomo_models/test_paper_ocp.py @@ -0,0 +1,318 @@ +# Copyright (C) 2026, SECQUOIA + +"""Tests for the paper-reference Pyomo OCP benchmark.""" + +import numpy as np +import pytest +from lyopronto.pyomo_models.paper_ocp import ( + PaperDiscretization, + PaperPrimaryDryingConfig, + classify_paper_policies, + create_paper_problem1_model, + derive_primary_drying_parameters, + generate_problem1_policy_initialization, + initialize_paper_problem1_from_trajectory, + interface_velocity, + load_upstream_matlab_trajectory, + product_resistance, + saturation_pressure, + solve_paper_problem1, + sublimation_flux, +) + +pyo = pytest.importorskip("pyomo.environ", reason="Pyomo not available") + +pytestmark = pytest.mark.pyomo + + +def _ipopt_available(): + try: + from idaes.core.solvers import get_solver + + return get_solver("ipopt").available() + except Exception: + try: + return pyo.SolverFactory("ipopt").available(False) + except Exception: + return False + + +def test_default_parameter_translation_matches_upstream_processing(): + config = PaperPrimaryDryingConfig() + derived = derive_primary_drying_parameters(config, n_z=20) + + assert np.isclose(derived.solution_density, 1018.86102386582) + assert np.isclose(derived.frozen_density, 936.7900511787848) + assert np.isclose(derived.frozen_heat_capacity, 2064.6) + assert np.isclose(derived.frozen_conductivity, 2.1438) + assert np.isclose(derived.frozen_diffusivity, 1.1084243905880022e-6) + assert np.isclose(derived.cross_section_area, 4.523893421169302e-4) + assert np.isclose(derived.product_height, 7.2124292981419705e-3) + assert len(derived.psi) == 20 + assert np.isclose(derived.dpsi, 1.0 / 19.0) + + +def test_equation_helpers_match_upstream_formula_values(): + config = PaperPrimaryDryingConfig() + derived = derive_primary_drying_parameters(config, n_z=20) + + assert np.isclose(saturation_pressure(228.0, config), 7.112217181578972) + assert np.isclose(saturation_pressure(243.0, config), 37.491783814828466) + + resistance = product_resistance(0.005, config) + expected_resistance = config.resistance_0 + config.resistance_1 * 0.005 / ( + 1.0 + config.resistance_2 * 0.005 + ) + assert np.isclose(resistance, expected_resistance) + + flux = sublimation_flux(228.0, 0.0, config) + expected_flux = (saturation_pressure(228.0, config) - 3.0) / config.resistance_0 + assert np.isclose(flux, expected_flux) + + velocity = interface_velocity(228.0, 0.0, config, derived) + assert np.isclose(velocity, flux / (derived.frozen_density - 215.0)) + + +def test_problem1_model_constructs_with_collocation(): + discretization = PaperDiscretization(n_z=5, nfe=4, ncp=2) + model = create_paper_problem1_model(discretization=discretization) + + assert len(list(model.z)) == 5 + assert len(list(model.t)) == 4 * 2 + 1 + assert hasattr(model, "T") + assert hasattr(model, "S") + assert hasattr(model, "Tb") + assert hasattr(model, "temperature_ode") + assert len(model.product_temperature_limit) == len(list(model.t)) + assert len(model.nonnegative_sublimation_flux) == len(list(model.t)) + assert hasattr(model, "terminal_drying") + assert hasattr(model, "objective") + + lb, ub = model.Tb[next(iter(model.t))].bounds + assert lb == 228.0 + assert ub == 273.0 + + +def test_problem1_model_initial_values_are_extractable(): + discretization = PaperDiscretization(n_z=5, nfe=3, ncp=2) + model = create_paper_problem1_model(discretization=discretization) + + t_points = sorted(model.t) + assert pyo.value(model.S[t_points[0]]) == 0.0 + assert pyo.value(model.S[t_points[-1]]) > 0.0 + assert pyo.value(model.t_final) == PaperPrimaryDryingConfig().problem1_time_guess + + +def test_problem1_model_constrains_sublimation_flux_nonnegative(): + config = PaperPrimaryDryingConfig() + discretization = PaperDiscretization(n_z=5, nfe=3, ncp=2) + model = create_paper_problem1_model(config=config, discretization=discretization) + + for t in model.t: + constraint = model.nonnegative_sublimation_flux[t] + assert constraint.lower == config.chamber_water_pressure + assert constraint.upper is None + + first_time = next(iter(model.t)) + model.T[0, first_time].set_value(discretization.temperature_lower_bound) + assert pyo.value(model.nonnegative_sublimation_flux[first_time].body) < ( + config.chamber_water_pressure + ) + + +def test_policy_classifier_detects_problem1_sequence(): + result = { + "states": { + "time_hr": np.array([0.0, 1.0, 2.0, 3.0, 4.0]), + "max_temperature_K": np.array([228.0, 235.0, 241.0, 242.8, 243.0]), + }, + "controls": { + "shelf_temperature_K": np.array([273.0, 273.0, 273.0, 265.0, 258.0]), + }, + "problem": { + "temperature_limit_K": 243.0, + "shelf_temperature_max_K": 273.0, + }, + } + + policies = classify_paper_policies(result, tolerances={"temperature_K": 0.3}) + + assert policies["segments"][0]["label"] == "policy_1_max_heat_input" + assert policies["segments"][1]["label"] == "policy_2_temperature_tracking" + assert policies["switch_times_hr"] == [3.0] + + +def test_policy_initialization_matches_upstream_policy1_event(): + trajectory = generate_problem1_policy_initialization( + discretization=PaperDiscretization(n_z=20), + n_time_points=80, + ) + + assert np.isclose( + trajectory["metrics"]["policy1_switch_time_hr"], + 2.363310733077, + atol=0.03, + ) + assert trajectory["policies"]["segments"][0]["label"] == "policy_1_max_heat_input" + assert trajectory["policies"]["segments"][1]["label"] == ( + "policy_2_temperature_tracking" + ) + assert trajectory["metrics"]["terminal_drying_fraction"] >= 0.994 + + +def test_initialize_model_from_policy_trajectory_sets_consistent_values(): + discretization = PaperDiscretization(n_z=5, nfe=4, ncp=2) + trajectory = generate_problem1_policy_initialization( + discretization=discretization, + n_time_points=80, + ) + model = create_paper_problem1_model(discretization=discretization) + + initialize_paper_problem1_from_trajectory(model, trajectory) + + t_points = sorted(model.t) + assert np.isclose( + pyo.value(model.t_final), + trajectory["metrics"]["drying_time_s"], + ) + assert np.isclose( + pyo.value(model.S[t_points[-1]]), + trajectory["states"]["interface_position_m"][-1], + atol=1e-7, + ) + assert pyo.value(model.Tb[t_points[0]]) == PaperPrimaryDryingConfig().shelf_temperature_max + assert pyo.value(model.Tb[t_points[-1]]) < PaperPrimaryDryingConfig().shelf_temperature_max + + +def test_load_upstream_matlab_trajectory_from_segment_file(tmp_path): + from scipy.io import savemat + + mat_path = tmp_path / "upstream_segment.mat" + t = np.array([0.0, 50.0, 100.0]) + y = np.array( + [ + [228.0, 228.1, 0.0], + [229.0, 230.0, 1.0e-4], + [230.0, 232.0, 2.0e-4], + ] + ) + savemat( + mat_path, + { + "t": t, + "y": y, + "Tb": np.full_like(t, 273.0), + "dSdt": np.array([1.0e-7, 1.1e-7, 1.2e-7]), + }, + ) + + trajectory = load_upstream_matlab_trajectory(mat_path) + + assert trajectory["metadata"]["source"] == "upstream_matlab" + assert trajectory["states"]["temperature_K"].shape == (3, 2) + assert np.allclose(trajectory["states"]["interface_position_m"], y[:, -1]) + assert np.allclose(trajectory["controls"]["shelf_temperature_K"], 273.0) + + +@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(): + discretization = PaperDiscretization( + n_z=5, + nfe=12, + ncp=3, + terminal_drying_fraction=0.995, + ) + result = solve_paper_problem1( + discretization=discretization, + solver_options={ + "max_iter": 2000, + "tol": 1.0e-5, + "acceptable_tol": 1.0e-4, + "print_level": 0, + }, + require_success=True, + ) + + metrics = result["metrics"] + policies = result["policies"] + labels = policies["labels"] + + assert metrics["terminal_gap_m"] <= 1.0e-7 + assert metrics["max_temperature_violation_K"] <= 1.0e-3 + assert metrics["shelf_lower_violation_K"] <= 1.0e-6 + assert metrics["shelf_upper_violation_K"] <= 1.0e-6 + assert np.isclose(metrics["drying_time_hr"], 6.19, atol=0.35) + assert "policy_1_max_heat_input" in labels + assert "policy_2_temperature_tracking" in labels + assert policies["segments"][0]["label"] == "policy_1_max_heat_input" + assert policies["segments"][1]["label"] == "policy_2_temperature_tracking" + assert np.isclose(policies["switch_times_hr"][0], 2.4, atol=0.35) + + +@pytest.mark.slow +@pytest.mark.skipif(not _ipopt_available(), reason="IPOPT solver not available") +def test_problem1_nz10_solve_matches_reference_policy_sequence(): + discretization = PaperDiscretization( + n_z=10, + nfe=12, + ncp=3, + terminal_drying_fraction=0.995, + ) + result = solve_paper_problem1( + discretization=discretization, + solver_options={ + "max_iter": 3000, + "tol": 1.0e-5, + "acceptable_tol": 1.0e-4, + "print_level": 0, + }, + require_success=True, + ) + + metrics = result["metrics"] + policies = result["policies"] + + assert result["metadata"]["status"] == "ok" + assert metrics["terminal_gap_m"] <= 1.0e-7 + assert metrics["max_temperature_violation_K"] <= 2.0e-6 + assert np.isclose(metrics["drying_time_hr"], 6.19, atol=0.08) + assert policies["segments"][0]["label"] == "policy_1_max_heat_input" + assert policies["segments"][1]["label"] == "policy_2_temperature_tracking" + assert np.isclose(policies["switch_times_hr"][0], 2.4, atol=0.12) + + +@pytest.mark.slow +@pytest.mark.skipif(not _ipopt_available(), reason="IPOPT solver not available") +def test_problem1_nz20_solve_matches_reference_policy_sequence(): + discretization = PaperDiscretization( + n_z=20, + nfe=12, + ncp=3, + terminal_drying_fraction=0.995, + ) + result = solve_paper_problem1( + discretization=discretization, + solver_options={ + "max_iter": 5000, + "max_cpu_time": 120, + "tol": 1.0e-5, + "acceptable_tol": 1.0e-3, + "acceptable_iter": 5, + "print_level": 0, + }, + require_success=True, + ) + + metrics = result["metrics"] + policies = result["policies"] + + assert result["metadata"]["status"] == "ok" + assert metrics["terminal_gap_m"] <= 1.0e-7 + assert metrics["max_temperature_violation_K"] <= 2.0e-6 + assert metrics["shelf_lower_violation_K"] <= 1.0e-6 + assert metrics["shelf_upper_violation_K"] <= 1.0e-6 + assert np.isclose(metrics["drying_time_hr"], 6.19, atol=0.08) + assert policies["segments"][0]["label"] == "policy_1_max_heat_input" + assert policies["segments"][1]["label"] == "policy_2_temperature_tracking" + assert np.isclose(policies["switch_times_hr"][0], 2.4, atol=0.12)