diff --git a/.github/workflows/pr-tests.yml b/.github/workflows/pr-tests.yml
index a8800fc..620ad9b 100644
--- a/.github/workflows/pr-tests.yml
+++ b/.github/workflows/pr-tests.yml
@@ -53,10 +53,10 @@ jobs:
       run: |
         if [ "${{ steps.mode.outputs.mode }}" == "fast" ]; then
           echo "Draft PR: Skipping notebook and pyomo tests for fast feedback"
-          pytest tests/ -n auto -v -m "not notebook and not pyomo"
+          pytest tests/ -n auto -v -m "not notebook and not pyomo" --ignore=tests/test_pyomo_models
         else
           echo "Ready PR: Running full scipy test suite (excluding notebook and pyomo)"
-          pytest tests/ -n auto -v -m "not notebook and not pyomo" --cov=lyopronto --cov-report=xml:coverage.xml
+          pytest tests/ -n auto -v -m "not notebook and not pyomo" --ignore=tests/test_pyomo_models --cov=lyopronto --cov-report=xml:coverage.xml
         fi
     
     - name: Upload coverage
diff --git a/.github/workflows/rundocs.yml b/.github/workflows/rundocs.yml
index e56fa60..c8fc41f 100644
--- a/.github/workflows/rundocs.yml
+++ b/.github/workflows/rundocs.yml
@@ -37,7 +37,7 @@ jobs:
         pip install -e . --no-build-isolation
     
     - name: Run tests (draft = fast, ready = coverage)
-      run: pytest tests/ -n auto -v -m "notebook" --cov=lyopronto --cov-report=xml --cov-report=term-missing
+      run: pytest tests/ -n auto -v -m "notebook" --ignore=tests/test_pyomo_models --cov=lyopronto --cov-report=xml --cov-report=term-missing
 # jobs:
 #   doctests:
 #     # env:
diff --git a/.github/workflows/slow-tests.yml b/.github/workflows/slow-tests.yml
index 20df2e0..11c4832 100644
--- a/.github/workflows/slow-tests.yml
+++ b/.github/workflows/slow-tests.yml
@@ -74,7 +74,7 @@ jobs:
           if [ "${{ inputs.include_pyomo }}" == "true" ]; then
             run_pytest pytest tests/ -n auto -v --cov=lyopronto --cov-report=xml --cov-report=term-missing
           else
-            run_pytest pytest tests/ -n auto -v -m "not pyomo" --cov=lyopronto --cov-report=xml --cov-report=term-missing
+            run_pytest pytest tests/ -n auto -v -m "not pyomo" --ignore=tests/test_pyomo_models --cov=lyopronto --cov-report=xml --cov-report=term-missing
           fi
         else
           echo "Running ONLY slow tests (marked with @pytest.mark.slow)"
@@ -82,7 +82,7 @@ jobs:
           if [ "${{ inputs.include_pyomo }}" == "true" ]; then
             run_pytest pytest tests/ -n auto -v -m "slow" --cov=lyopronto --cov-report=xml --cov-report=term-missing
           else
-            run_pytest pytest tests/ -n auto -v -m "slow and not pyomo" --cov=lyopronto --cov-report=xml --cov-report=term-missing
+            run_pytest pytest tests/ -n auto -v -m "slow and not pyomo" --ignore=tests/test_pyomo_models --cov=lyopronto --cov-report=xml --cov-report=term-missing
           fi
         fi
     
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index e3be451..04b8359 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -37,7 +37,7 @@ jobs:
     - name: Run SciPy tests (excluding Pyomo tests)
       run: |
         echo "Running SciPy test suite (excluding Pyomo tests)"
-        pytest tests/ -n auto -v -m "not pyomo" --cov=lyopronto --cov-report=xml --cov-report=term-missing
+        pytest tests/ -n auto -v -m "not pyomo" --ignore=tests/test_pyomo_models --cov=lyopronto --cov-report=xml --cov-report=term-missing
     
     - name: Upload coverage to Codecov
       uses: codecov/codecov-action@v4
diff --git a/lyopronto/pyomo_models/single_step.py b/lyopronto/pyomo_models/single_step.py
new file mode 100644
index 0000000..d52d677
--- /dev/null
+++ b/lyopronto/pyomo_models/single_step.py
@@ -0,0 +1,427 @@
+"""Single time-step Pyomo optimization for lyophilization primary drying.
+
+This module provides a Pyomo-based NLP formulation that replicates one time step
+of the scipy sequential optimization approach. It solves for optimal chamber
+pressure (Pch) and shelf temperature (Tsh) at a given dried cake length (Lck).
+
+Physics Corrections (Jan 2025):
+    - Energy balance: Corrected to match multi-period model (Q_shelf = Q_sublimation)
+    - Vial bottom temp: Added proper constraint for frozen layer temperature gradient
+    - Matches corrected multi-period formulation
+
+The model includes:
+    - 7 decision variables: Pch, Tsh, Tsub, Tbot, Psub, dmdt, Kv
+    - 5 equality constraints: vapor pressure, sublimation rate, energy balance,
+      vial bottom temperature, vial heat transfer
+    - 2 inequality constraints: equipment capability, product temperature limit
+    - Objective: maximize sublimation driving force (minimize Pch - Psub)
+"""
+
+# LyoPRONTO, a vial-scale lyophilization process simulator
+# Nonlinear optimization
+# Copyright (C) 2025, David E. Bernal Neira
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import pyomo.environ as pyo
+
+from .. import constant
+
+
+def create_single_step_model(
+    vial,
+    product,
+    ht,
+    Lpr0,
+    Lck,
+    Pch_bounds=(0.05, 0.5),
+    Tsh_bounds=(-50, 50),
+    eq_cap=None,
+    nVial=None,
+    apply_scaling=True,
+):
+    """Create Pyomo model for single time-step lyophilization optimization.
+
+    This function constructs a ConcreteModel that represents the physics and
+    constraints of lyophilization at a single time step. The model finds optimal
+    chamber pressure and shelf temperature to maximize sublimation rate while
+    respecting equipment and product temperature constraints.
+
+    Args:
+        vial (dict): Vial geometry with keys:
+            - 'Av' (float): Vial area [cm²]
+            - 'Ap' (float): Product area [cm²]
+        product (dict): Product properties with keys:
+            - 'R0' (float): Base product resistance [cm²·hr·Torr/g]
+            - 'A1' (float): Resistance parameter [cm·hr·Torr/g]
+            - 'A2' (float): Resistance parameter [1/cm]
+            - 'T_pr_crit' (float): Critical product temperature [°C]
+        ht (dict): Heat transfer parameters with keys:
+            - 'KC' (float): Vial heat transfer parameter [cal/s/K/cm²]
+            - 'KP' (float): Vial heat transfer parameter [cal/s/K/cm²/Torr]
+            - 'KD' (float): Vial heat transfer parameter [1/Torr]
+        Lpr0 (float): Initial product length [cm]
+        Lck (float): Current dried cake length [cm]
+        Pch_bounds (tuple, optional): (min, max) for chamber pressure [Torr].
+            Default: (0.05, 0.5)
+        Tsh_bounds (tuple, optional): (min, max) for shelf temperature [°C].
+            Default: (-50, 50)
+        eq_cap (dict, optional): Equipment capability parameters with keys:
+            - 'a' (float): Equipment capability parameter [kg/hr]
+            - 'b' (float): Equipment capability parameter [kg/hr/Torr]
+            If None, equipment constraint is not enforced.
+        nVial (int, optional): Number of vials in batch. Required if eq_cap provided.
+
+    Returns:
+        pyo.ConcreteModel: Pyomo model ready to solve with variables, constraints,
+            and objective defined.
+
+    Notes:
+        The model uses the following physics equations (corrected Jan 2025):
+        - Vapor pressure: Psub = 2.698e10 * exp(-6144.96/(Tsub + 273.15))
+        - Product resistance: Rp = R0 + A1*Lck/(1 + A2*Lck)
+        - Vial heat transfer: Kv = KC + KP*Pch/(1 + KD*Pch)
+        - Sublimation rate: dmdt = Ap/Rp * (Psub - Pch) / kg_To_g
+        - Energy balance: Kv*Av*(Tsh - Tbot) = dHs*dmdt*kg_To_g/hr_To_s
+        - Vial bottom temp: Tbot = Tsub + (Lpr0-Lck)*dHs*dmdt*kg_To_g/hr_To_s/(Ap*k_ice)
+
+    Examples:
+        >>> from lyopronto import functions
+        >>> vial = {'Av': 3.80, 'Ap': 3.14}
+        >>> product = {'R0': 1.4, 'A1': 16.0, 'A2': 0.0, 'T_pr_crit': -5.0}
+        >>> ht = {'KC': 2.75e-4, 'KP': 8.93e-4, 'KD': 0.46}
+        >>> Lpr0 = functions.Lpr0_FUN(2.0, 3.14, 0.05)
+        >>> Lck = 0.5  # Half dried
+        >>> model = create_single_step_model(vial, product, ht, Lpr0, Lck)
+        >>> # Now solve with solve_single_step(model)
+    """
+    model = pyo.ConcreteModel()
+
+    # ==================== Parameters (Fixed for this step) ====================
+    model.Lpr0 = pyo.Param(initialize=Lpr0)
+    model.Lck = pyo.Param(initialize=Lck)
+    model.Av = pyo.Param(initialize=vial["Av"])
+    model.Ap = pyo.Param(initialize=vial["Ap"])
+    model.R0 = pyo.Param(initialize=product["R0"])
+    model.A1 = pyo.Param(initialize=product["A1"])
+    model.A2 = pyo.Param(initialize=product["A2"])
+    model.T_crit = pyo.Param(initialize=product["T_pr_crit"])
+    model.KC = pyo.Param(initialize=ht["KC"])
+    model.KP = pyo.Param(initialize=ht["KP"])
+    model.KD = pyo.Param(initialize=ht["KD"])
+
+    # Physical constants
+    model.kg_To_g = pyo.Param(initialize=constant.kg_To_g)
+    model.hr_To_s = pyo.Param(initialize=constant.hr_To_s)
+    model.k_ice = pyo.Param(initialize=constant.k_ice)
+    model.dHs = pyo.Param(initialize=constant.dHs)
+
+    # ==================== Decision Variables ====================
+    # Chamber pressure [Torr]
+    model.Pch = pyo.Var(domain=pyo.NonNegativeReals, bounds=Pch_bounds)
+
+    # Shelf temperature [°C]
+    model.Tsh = pyo.Var(domain=pyo.Reals, bounds=Tsh_bounds)
+
+    # Sublimation front temperature [°C] - always below freezing
+    model.Tsub = pyo.Var(domain=pyo.Reals, bounds=(-60, 0))
+
+    # Vial bottom temperature [°C]
+    model.Tbot = pyo.Var(domain=pyo.Reals, bounds=(-60, 50))
+
+    # Vapor pressure at sublimation front [Torr]
+    model.Psub = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1e-6, 10))
+
+    # Log of vapor pressure (for numerical stability in exponential constraint)
+    model.log_Psub = pyo.Var(
+        domain=pyo.Reals, bounds=(-14, 2.5)
+    )  # log(1e-6) to log(10)
+
+    # Sublimation rate [kg/hr] - must be non-negative
+    model.dmdt = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0, 10))
+
+    # Vial heat transfer coefficient [cal/s/K/cm²]
+    model.Kv = pyo.Var(domain=pyo.PositiveReals, bounds=(1e-6, 1e-2))
+
+    # ==================== Derived Quantities (Expressions) ====================
+    # Product resistance [cm²·hr·Torr/g]
+    model.Rp = pyo.Expression(
+        expr=model.R0 + model.A1 * model.Lck / (1 + model.A2 * model.Lck)
+    )
+
+    # ==================== Equality Constraints ====================
+    # C1: Vapor pressure at sublimation front (Antoine equation)
+    # Using log transformation for numerical stability:
+    # log(Psub) = log(2.698e10) - 6144.96/(Tsub + 273.15)
+    model.vapor_pressure_log = pyo.Constraint(
+        expr=model.log_Psub == pyo.log(2.698e10) - 6144.96 / (model.Tsub + 273.15)
+    )
+
+    # C1b: Link log_Psub to Psub
+    model.vapor_pressure_exp = pyo.Constraint(
+        expr=model.Psub == pyo.exp(model.log_Psub)
+    )
+
+    # C2: Sublimation rate from mass transfer
+    # dmdt = Ap/Rp * (Psub - Pch) / kg_To_g
+    model.sublimation_rate = pyo.Constraint(
+        expr=model.dmdt
+        == model.Ap / model.Rp / model.kg_To_g * (model.Psub - model.Pch)
+    )
+
+    # C3: Energy balance - heat from shelf equals heat for sublimation
+    # Q_shelf = Kv * Av * (Tsh - Tbot)
+    # Q_sub = dHs * dmdt * kg_To_g / hr_To_s
+    # This matches the corrected multi-period energy_balance constraint
+    model.energy_balance = pyo.Constraint(
+        expr=model.Kv * model.Av * (model.Tsh - model.Tbot)
+        == model.dHs * model.dmdt * model.kg_To_g / model.hr_To_s
+    )
+
+    # C4: Vial bottom temperature from conduction through frozen layer
+    # Tbot = Tsub + ΔT_frozen
+    # where ΔT_frozen = (Lpr0 - Lck) * Q_sub / (Ap * k_ice)
+    # This matches the corrected multi-period vial_bottom_temp constraint
+    model.vial_bottom_temp = pyo.Constraint(
+        expr=model.Tbot
+        == model.Tsub
+        + (model.Lpr0 - model.Lck)
+        * model.dHs
+        * model.dmdt
+        * model.kg_To_g
+        / model.hr_To_s
+        / (model.Ap * model.k_ice)
+    )
+
+    # C5: Vial heat transfer coefficient
+    # Kv = KC + KP * Pch / (1 + KD * Pch)
+    model.kv_calc = pyo.Constraint(
+        expr=model.Kv == model.KC + model.KP * model.Pch / (1 + model.KD * model.Pch)
+    )
+
+    # ==================== Inequality Constraints ====================
+    # Product temperature limit: Tbot ≤ T_crit
+    model.temp_limit = pyo.Constraint(expr=model.Tbot <= model.T_crit)
+
+    # Equipment capability limit (optional)
+    if eq_cap is not None and nVial is not None:
+        model.a_eq = pyo.Param(initialize=eq_cap["a"])
+        model.b_eq = pyo.Param(initialize=eq_cap["b"])
+        model.nVial = pyo.Param(initialize=nVial)
+
+        # a + b*Pch - nVial*dmdt ≥ 0
+        model.equipment_capability = pyo.Constraint(
+            expr=model.a_eq + model.b_eq * model.Pch - model.nVial * model.dmdt >= 0
+        )
+
+    # ==================== Objective Function ====================
+    # Minimize (Pch - Psub) to maximize sublimation driving force (Psub - Pch)
+    model.obj = pyo.Objective(expr=model.Pch - model.Psub, sense=pyo.minimize)
+
+    # ==================== Scaling (Optional) ====================
+    if apply_scaling:
+        from . import utils
+
+        scaling_factors = {
+            "Tsub": 0.01,
+            "Tbot": 0.01,
+            "Tsh": 0.01,
+            "Pch": 10,
+            "Psub": 10,
+            "log_Psub": 1.0,  # log values are already scaled
+            "Kv": 1e4,
+            "dmdt": 1.0,
+        }
+        utils.add_scaling_suffix(model, scaling_factors)
+
+    return model
+
+
+def solve_single_step(model, solver="ipopt", tee=False, warmstart_data=None):
+    """Solve Pyomo single-step model and extract results.
+
+    Args:
+        model (pyo.ConcreteModel): Pyomo model created by create_single_step_model()
+        solver (str, optional): Solver name. Default: 'ipopt'
+        tee (bool, optional): If True, print solver output. Default: False
+        warmstart_data (dict, optional): Initial values for variables with keys:
+            'Pch', 'Tsh', 'Tsub', 'Tbot', 'Psub', 'dmdt', 'Kv'.
+            If provided, these values are used to initialize the model.
+
+    Returns:
+        dict: Solution dictionary with keys:
+            - 'status' (str): Solver termination condition
+            - 'Pch' (float): Chamber pressure [Torr]
+            - 'Tsh' (float): Shelf temperature [°C]
+            - 'Tsub' (float): Sublimation front temperature [°C]
+            - 'Tbot' (float): Vial bottom temperature [°C]
+            - 'Psub' (float): Vapor pressure [Torr]
+            - 'dmdt' (float): Sublimation rate [kg/hr]
+            - 'Kv' (float): Vial heat transfer coefficient [cal/s/K/cm²]
+            - 'Rp' (float): Product resistance [cm²·hr·Torr/g]
+            - 'obj' (float): Objective value
+
+    Raises:
+        RuntimeError: If solver is not available or fails to solve.
+
+    Notes:
+        For IPOPT solver, the following options are used:
+        - max_iter: 3000
+        - tol: 1e-6
+        - mu_strategy: 'adaptive'
+        - print_level: 5 (if tee=True), 0 (if tee=False)
+
+        If IPOPT is not in PATH, will attempt to use IDAES-provided IPOPT.
+
+    Examples:
+        >>> model = create_single_step_model(vial, product, ht, Lpr0, Lck)
+        >>> solution = solve_single_step(model, tee=True)
+        >>> print(f"Optimal Pch: {solution['Pch']:.4f} Torr")
+        >>> print(f"Optimal Tsh: {solution['Tsh']:.2f} °C")
+    """
+    # Apply warmstart values if provided
+    if warmstart_data is not None:
+        if "Pch" in warmstart_data:
+            model.Pch.set_value(warmstart_data["Pch"])
+        if "Tsh" in warmstart_data:
+            model.Tsh.set_value(warmstart_data["Tsh"])
+        if "Tsub" in warmstart_data:
+            model.Tsub.set_value(warmstart_data["Tsub"])
+        if "Tbot" in warmstart_data:
+            model.Tbot.set_value(warmstart_data["Tbot"])
+        if "Psub" in warmstart_data:
+            model.Psub.set_value(warmstart_data["Psub"])
+        if "log_Psub" in warmstart_data:
+            model.log_Psub.set_value(warmstart_data["log_Psub"])
+        if "dmdt" in warmstart_data:
+            model.dmdt.set_value(warmstart_data["dmdt"])
+        if "Kv" in warmstart_data:
+            model.Kv.set_value(warmstart_data["Kv"])
+
+    # Create solver - try IDAES first if available, then standard Pyomo
+    if solver.lower() == "ipopt":
+        try:
+            from idaes.core.solvers import get_solver
+
+            opt = get_solver("ipopt")
+        except (ImportError, Exception):
+            # Fall back to standard Pyomo solver
+            opt = pyo.SolverFactory(solver)
+    else:
+        opt = pyo.SolverFactory(solver)
+
+    # Configure IPOPT options
+    if solver.lower() == "ipopt":
+        opt.options["max_iter"] = 3000
+        opt.options["tol"] = 1e-6
+        opt.options["mu_strategy"] = "adaptive"
+        opt.options["print_level"] = 5 if tee else 0
+
+    # Solve
+    results = opt.solve(model, tee=tee)
+
+    # Check convergence
+    termination = results.solver.termination_condition
+    status_str = str(termination)
+
+    if termination not in [
+        pyo.TerminationCondition.optimal,
+        pyo.TerminationCondition.locallyOptimal,
+    ]:
+        print(f"WARNING: Solver status: {termination}")
+
+    # Extract solution
+    solution = {
+        "status": status_str,
+        "Pch": pyo.value(model.Pch),
+        "Tsh": pyo.value(model.Tsh),
+        "Tsub": pyo.value(model.Tsub),
+        "Tbot": pyo.value(model.Tbot),
+        "Psub": pyo.value(model.Psub),
+        "log_Psub": pyo.value(model.log_Psub),
+        "dmdt": pyo.value(model.dmdt),
+        "Kv": pyo.value(model.Kv),
+        "Rp": pyo.value(model.Rp),
+        "obj": pyo.value(model.obj),
+    }
+
+    return solution
+
+
+def optimize_single_step(
+    vial,
+    product,
+    ht,
+    Lpr0,
+    Lck,
+    Pch_bounds=(0.05, 0.5),
+    Tsh_bounds=(-50, 50),
+    eq_cap=None,
+    nVial=None,
+    warmstart_data=None,
+    solver="ipopt",
+    tee=False,
+    apply_scaling=True,
+):
+    """Convenience function to create and solve single-step model in one call.
+
+    This function combines create_single_step_model() and solve_single_step()
+    for ease of use when you don't need to inspect the model structure.
+
+    Args:
+        vial (dict): Vial geometry (see create_single_step_model)
+        product (dict): Product properties (see create_single_step_model)
+        ht (dict): Heat transfer parameters (see create_single_step_model)
+        Lpr0 (float): Initial product length [cm]
+        Lck (float): Current dried cake length [cm]
+        Pch_bounds (tuple): Chamber pressure bounds [Torr]
+        Tsh_bounds (tuple): Shelf temperature bounds [°C]
+        eq_cap (dict): Equipment capability parameters
+        nVial (int): Number of vials
+        warmstart_data (dict): Initial variable values
+        solver (str): Solver name (default: 'ipopt')
+        tee (bool): Print solver output (default: False)
+        apply_scaling (bool): Apply variable scaling (default: True)
+
+    Returns:
+        dict: Solution dictionary (see solve_single_step)
+
+    Examples:
+        >>> solution = optimize_single_step(
+        ...     vial={'Av': 3.8, 'Ap': 3.14},
+        ...     product={'R0': 1.4, 'A1': 16.0, 'A2': 0.0, 'T_pr_crit': -5.0},
+        ...     ht={'KC': 2.75e-4, 'KP': 8.93e-4, 'KD': 0.46},
+        ...     Lpr0=1.5,
+        ...     Lck=0.5,
+        ...     tee=True
+        ... )
+    """
+    # Create and solve
+    model = create_single_step_model(
+        vial,
+        product,
+        ht,
+        Lpr0,
+        Lck,
+        Pch_bounds=Pch_bounds,
+        Tsh_bounds=Tsh_bounds,
+        eq_cap=eq_cap,
+        nVial=nVial,
+        apply_scaling=apply_scaling,
+    )
+    solution = solve_single_step(
+        model, solver=solver, tee=tee, warmstart_data=warmstart_data
+    )
+
+    return solution
diff --git a/lyopronto/pyomo_models/utils.py b/lyopronto/pyomo_models/utils.py
new file mode 100644
index 0000000..7f4e94f
--- /dev/null
+++ b/lyopronto/pyomo_models/utils.py
@@ -0,0 +1,254 @@
+"""Utility functions for Pyomo model initialization, scaling, and result extraction.
+
+This module provides helper functions to bridge scipy and Pyomo implementations,
+including warmstarting Pyomo models from scipy solutions and converting results
+to standard output formats.
+"""
+
+# LyoPRONTO, a vial-scale lyophilization process simulator
+# Nonlinear optimization
+# Copyright (C) 2025, David E. Bernal Neira
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import numpy as np
+
+from .. import constant, functions
+
+
+def initialize_from_scipy(scipy_output, time_index, vial, product, Lpr0, ht=None):
+    """Create warmstart data dictionary from scipy optimization output.
+
+    This function extracts values from a scipy optimization result array
+    and formats them as a dictionary suitable for initializing Pyomo variables.
+
+    Args:
+        scipy_output (np.ndarray): Output from opt_Pch_Tsh.dry() or similar,
+            with shape (n_steps, 7) and columns:
+            [time, Tsub, Tbot, Tsh, Pch_mTorr, flux, percent_dried (0-100)]
+        time_index (int): Index of the time step to extract (0-based)
+        vial (dict): Vial geometry with 'Av', 'Ap' keys
+        product (dict): Product properties with 'R0', 'A1', 'A2' keys
+        Lpr0 (float): Initial product length [cm]
+        ht (dict, optional): Heat transfer parameters with 'KC', 'KP', 'KD' keys.
+            If provided, Kv will be computed accurately.
+
+    Returns:
+        dict: Warmstart data with keys: 'Pch', 'Tsh', 'Tsub', 'Tbot', 'Psub',
+            'log_Psub', 'dmdt', 'Kv'
+
+    Notes:
+        - Pch is converted from mTorr to Torr (divides by 1000)
+        - Lck is calculated from percent_dried (column 6 is 0-100%, converted to fraction)
+        - Derived quantities (Psub, dmdt, Kv) are computed from physics functions
+
+    Examples:
+        >>> from lyopronto import opt_Pch_Tsh
+        >>> scipy_out = opt_Pch_Tsh.dry(vial, product, ht, Pch, Tsh, dt, eq_cap, nVial)
+        >>> warmstart = initialize_from_scipy(scipy_out, 10, vial, product, Lpr0)
+        >>> # Use warmstart dict to initialize Pyomo model
+    """
+    # Extract values from scipy output
+    # Columns: [time, Tsub, Tbot, Tsh, Pch_mTorr, flux, percent_dried]
+    Tsub = scipy_output[time_index, 1]
+    Tbot = scipy_output[time_index, 2]
+    Tsh = scipy_output[time_index, 3]
+    Pch = scipy_output[time_index, 4] / constant.Torr_to_mTorr  # mTorr → Torr
+    percent_dried = scipy_output[time_index, 6]  # 0-100 percentage
+
+    # Calculate derived quantities
+    frac_dried = percent_dried / 100.0  # Convert percentage to fraction (0-1)
+    Lck = frac_dried * Lpr0  # Current dried cake length [cm]
+    Rp = functions.Rp_FUN(Lck, product["R0"], product["A1"], product["A2"])
+    Psub = functions.Vapor_pressure(Tsub)
+    dmdt = functions.sub_rate(vial["Ap"], Rp, Tsub, Pch)
+
+    # Calculate Kv from heat transfer parameters if available
+    if ht is not None:
+        Kv = functions.Kv_FUN(ht["KC"], ht["KP"], ht["KD"], Pch)
+    else:
+        # Use typical value as fallback
+        Kv = 5e-4  # Typical value [cal/s/K/cm²]
+
+    warmstart_data = {
+        "Pch": Pch,
+        "Tsh": Tsh,
+        "Tsub": Tsub,
+        "Tbot": Tbot,
+        "Psub": Psub,
+        "log_Psub": np.log(max(Psub, 1e-10)),  # Add log for stability
+        "dmdt": max(0.0, dmdt),  # Ensure non-negative
+        "Kv": Kv,
+    }
+
+    return warmstart_data
+
+
+def extract_solution_to_array(solution, time):
+    """Convert Pyomo solution dict to standard output array format.
+
+    This function formats a Pyomo solution to match the scipy output format
+    for consistency and comparison.
+
+    Args:
+        solution (dict): Solution from solve_single_step() with keys:
+            'Pch', 'Tsh', 'Tsub', 'Tbot', 'dmdt', etc.
+        time (float): Time value for this step [hr]
+
+    Returns:
+        np.ndarray: Array of shape (7,) with columns:
+            [time, Tsub, Tbot, Tsh, Pch_mTorr, flux, percent_dried]
+
+    Notes:
+        - Pch is converted from Torr to mTorr
+        - flux is dmdt normalized by product area
+        - percent_dried must be computed externally (requires Lck and Lpr0)
+
+    Examples:
+        >>> solution = solve_single_step(model)
+        >>> output_row = extract_solution_to_array(solution, time=0.5)
+    """
+    # Note: This is a simplified version
+    # In full implementation, would need vial['Ap'] and Lck/Lpr0 for complete conversion
+    output = np.array(
+        [
+            time,
+            solution["Tsub"],
+            solution["Tbot"],
+            solution["Tsh"],
+            solution["Pch"] * constant.Torr_to_mTorr,  # Torr → mTorr
+            solution["dmdt"],  # Note: needs conversion to flux [kg/hr/m²]
+            0.0,  # percent_dried - needs to be computed externally
+        ]
+    )
+
+    return output
+
+
+def add_scaling_suffix(model, variable_scales=None):
+    """Add scaling factors to Pyomo model for improved numerical conditioning.
+
+    Scaling can significantly improve solver convergence by ensuring all
+    variables and constraints have similar magnitudes.
+
+    Args:
+        model (pyo.ConcreteModel): Pyomo model to add scaling to
+        variable_scales (dict, optional): Custom scaling factors. If None,
+            uses default scales based on typical variable magnitudes.
+            Keys are variable names ('Tsub', 'Pch', etc.), values are
+            scaling factors.
+
+    Returns:
+        None: Modifies model in place by adding scaling_factor Suffix
+
+    Notes:
+        Default scaling factors:
+        - Temperature variables (Tsub, Tbot, Tsh): 0.01 (typical ~-20°C)
+        - Pressure variables (Pch, Psub): 10 (typical ~0.1 Torr)
+        - Kv: 1e4 (typical ~1e-4)
+        - dmdt: 1.0
+
+    Examples:
+        >>> model = create_single_step_model(...)
+        >>> add_scaling_suffix(model)  # Use defaults
+        >>> # Or with custom scales:
+        >>> add_scaling_suffix(model, {'Pch': 5, 'Tsh': 0.02})
+    """
+    import pyomo.environ as pyo
+
+    # Default scales
+    default_scales = {
+        "Tsub": 0.01,  # Typical value ~-20°C
+        "Tbot": 0.01,
+        "Tsh": 0.01,
+        "Pch": 10,  # Typical value ~0.1 Torr
+        "Psub": 10,
+        "Kv": 1e4,  # Typical value ~1e-4
+        "dmdt": 1.0,
+    }
+
+    # Use custom scales if provided, otherwise defaults
+    scales = variable_scales if variable_scales is not None else default_scales
+
+    # Create scaling suffix
+    model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT)
+
+    # Apply scaling factors
+    for var_name, scale in scales.items():
+        if hasattr(model, var_name):
+            var = getattr(model, var_name)
+            model.scaling_factor[var] = scale
+
+
+def check_solution_validity(solution, tol=1e-3):
+    """Validate that a Pyomo solution satisfies physical constraints.
+
+    Args:
+        solution (dict): Solution dictionary from solve_single_step()
+        tol (float, optional): Tolerance for constraint violations. Default: 1e-3
+
+    Returns:
+        tuple: (is_valid, violations) where:
+            - is_valid (bool): True if all checks pass
+            - violations (list): List of violation messages
+
+    Notes:
+        Checks performed:
+        - Temperature ordering: Tsub ≤ Tbot ≤ Tsh
+        - Sublimation temperature below freezing: Tsub < 0
+        - Positive pressures and rates
+        - Vapor pressure consistency
+
+    Examples:
+        >>> solution = solve_single_step(model)
+        >>> is_valid, violations = check_solution_validity(solution)
+        >>> if not is_valid:
+        ...     print("Violations:", violations)
+    """
+    violations = []
+
+    # Temperature ordering
+    if solution["Tsub"] > solution["Tbot"] + tol:
+        violations.append(
+            f"Tsub ({solution['Tsub']:.2f}) > Tbot ({solution['Tbot']:.2f})"
+        )
+
+    if solution["Tbot"] > solution["Tsh"] + tol:
+        violations.append(
+            f"Tbot ({solution['Tbot']:.2f}) > Tsh ({solution['Tsh']:.2f})"
+        )
+
+    # Sublimation temperature below freezing
+    if solution["Tsub"] > 0 + tol:
+        violations.append(f"Tsub ({solution['Tsub']:.2f}) above freezing")
+
+    # Positive values
+    if solution["Pch"] < -tol:
+        violations.append(f"Negative Pch: {solution['Pch']:.4f}")
+
+    if solution["Psub"] < -tol:
+        violations.append(f"Negative Psub: {solution['Psub']:.4f}")
+
+    if solution["dmdt"] < -tol:
+        violations.append(f"Negative dmdt: {solution['dmdt']:.6f}")
+
+    # Driving force check (Psub should be > Pch for sublimation)
+    if solution["Psub"] < solution["Pch"] - tol:
+        violations.append(
+            f"Psub ({solution['Psub']:.4f}) < Pch ({solution['Pch']:.4f})"
+        )
+
+    is_valid = len(violations) == 0
+
+    return is_valid, violations
diff --git a/tests/test_pyomo_models/README.md b/tests/test_pyomo_models/README.md
new file mode 100644
index 0000000..7bb5618
--- /dev/null
+++ b/tests/test_pyomo_models/README.md
@@ -0,0 +1,29 @@
+# Pyomo Models Test Suite
+
+Tests for the Pyomo-based optimization models in `lyopronto/pyomo_models/`.
+
+## Test Organization
+
+| Category | Files | Purpose |
+|----------|-------|---------|
+| **Model Tests** | `test_model_*.py` | Model creation, structure, warmstart |
+| **Optimizer Tests** | `test_optimizer_*.py` | Optimizer functions (Tsh, Pch, Pch_Tsh) |
+| **Infrastructure** | `test_parameter_validation.py`, `test_warmstart.py`, `test_staged_solve.py` | Validation, warmstart, staged solve |
+
+## Running Tests
+
+```bash
+# Run all Pyomo tests
+pytest tests/test_pyomo_models/ -v
+
+# Run with coverage
+pytest tests/test_pyomo_models/ --cov=lyopronto.pyomo_models
+
+# Skip if Pyomo not installed
+pytest tests/test_pyomo_models/ -v  # Tests auto-skip if Pyomo unavailable
+```
+
+## Known Limitations
+
+- `test_block_triangularization` in `test_model_multi_period.py` is xfailed due to DOF from initial conditions (structural analysis limitation, not a functional bug)
+- Some tests require optional dependencies (PyNumero, incidence_analysis)
diff --git a/tests/test_pyomo_models/test_model_single_step.py b/tests/test_pyomo_models/test_model_single_step.py
new file mode 100644
index 0000000..6ac8fe0
--- /dev/null
+++ b/tests/test_pyomo_models/test_model_single_step.py
@@ -0,0 +1,307 @@
+"""Tests for Pyomo single time-step model.
+
+This module tests the single time-step optimization model (from model.single_step)
+against the scipy baseline to ensure correctness and consistency.
+"""
+
+# LyoPRONTO, a vial-scale lyophilization process simulator
+# Nonlinear optimization
+# Copyright (C) 2025, David E. Bernal Neira
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import numpy as np
+import pytest
+from lyopronto import functions
+from lyopronto.pyomo_models import single_step, utils
+
+# Try to import pyomo - skip tests if not available
+try:
+    import pyomo.environ as pyo
+
+    PYOMO_AVAILABLE = True
+except ImportError:
+    PYOMO_AVAILABLE = False
+
+pytestmark = [
+    pytest.mark.pyomo,
+    pytest.mark.skipif(not PYOMO_AVAILABLE, reason="Pyomo not installed"),
+]
+
+
+class TestSingleStepModel:
+    """Tests for single time-step model creation and structure."""
+
+    def test_model_creation_basic(self, standard_vial, standard_product, standard_ht):
+        """Test that model can be created with standard inputs."""
+        Lpr0 = functions.Lpr0_FUN(2.0, standard_vial["Ap"], standard_product["cSolid"])
+        Lck = 0.5  # Half dried
+
+        model = single_step.create_single_step_model(
+            standard_vial, standard_product, standard_ht, Lpr0, Lck
+        )
+
+        # Check model type
+        assert isinstance(model, pyo.ConcreteModel)
+
+        # Check key variables exist
+        assert hasattr(model, "Pch")
+        assert hasattr(model, "Tsh")
+        assert hasattr(model, "Tsub")
+        assert hasattr(model, "Tbot")
+        assert hasattr(model, "Psub")
+        assert hasattr(model, "dmdt")
+        assert hasattr(model, "Kv")
+
+        # Check constraints exist (updated Jan 2025 for corrected physics)
+        assert hasattr(model, "vapor_pressure_log")
+        assert hasattr(model, "vapor_pressure_exp")
+        assert hasattr(model, "sublimation_rate")
+        assert hasattr(model, "energy_balance")  # Renamed from heat_balance
+        assert hasattr(model, "vial_bottom_temp")  # Renamed from shelf_temp
+        assert hasattr(model, "kv_calc")
+        assert hasattr(model, "temp_limit")
+
+        # Check log_Psub variable exists
+        assert hasattr(model, "log_Psub")
+
+        # Check objective exists
+        assert hasattr(model, "obj")
+
+    def test_model_with_equipment_constraint(
+        self, standard_vial, standard_product, standard_ht
+    ):
+        """Test model creation with equipment capability constraint."""
+        Lpr0 = functions.Lpr0_FUN(2.0, standard_vial["Ap"], standard_product["cSolid"])
+        Lck = 0.5
+        eq_cap = {"a": -0.182, "b": 11.7}
+        nVial = 398
+
+        model = single_step.create_single_step_model(
+            standard_vial,
+            standard_product,
+            standard_ht,
+            Lpr0,
+            Lck,
+            eq_cap=eq_cap,
+            nVial=nVial,
+        )
+
+        # Check equipment constraint exists
+        assert hasattr(model, "equipment_capability")
+        assert hasattr(model, "a_eq")
+        assert hasattr(model, "b_eq")
+        assert hasattr(model, "nVial")
+
+    def test_variable_bounds(self, standard_vial, standard_product, standard_ht):
+        """Test that variable bounds are correctly set."""
+        Lpr0 = functions.Lpr0_FUN(2.0, standard_vial["Ap"], standard_product["cSolid"])
+        Lck = 0.5
+        Pch_bounds = (0.1, 0.3)
+        Tsh_bounds = (-40, 30)
+
+        model = single_step.create_single_step_model(
+            standard_vial,
+            standard_product,
+            standard_ht,
+            Lpr0,
+            Lck,
+            Pch_bounds=Pch_bounds,
+            Tsh_bounds=Tsh_bounds,
+        )
+
+        # Check bounds
+        assert model.Pch.bounds == Pch_bounds
+        assert model.Tsh.bounds == Tsh_bounds
+        assert model.Tsub.bounds == (-60, 0)
+        assert model.dmdt.bounds[0] == 0  # Lower bound must be non-negative
+
+
+class TestSingleStepSolver:
+    """Tests for solving Pyomo single-step model."""
+
+    @pytest.mark.slow
+    def test_solve_basic(self, standard_vial, standard_product, standard_ht):
+        """Test that model can be solved successfully."""
+        Lpr0 = functions.Lpr0_FUN(2.0, standard_vial["Ap"], standard_product["cSolid"])
+        Lck = 0.5
+
+        model = single_step.create_single_step_model(
+            standard_vial, standard_product, standard_ht, Lpr0, Lck
+        )
+
+        solution = single_step.solve_single_step(model, tee=False)
+
+        # Check solution structure
+        assert "status" in solution
+        assert "Pch" in solution
+        assert "Tsh" in solution
+        assert "Tsub" in solution
+        assert "Tbot" in solution
+        assert "Psub" in solution
+        assert "dmdt" in solution
+        assert "Kv" in solution
+
+        # Check physical validity
+        assert solution["Tsub"] < 0, "Sublimation temp should be below freezing"
+        assert solution["Tsub"] <= solution["Tbot"], "Tsub should be <= Tbot"
+        assert solution["Tbot"] <= solution["Tsh"], "Tbot should be <= Tsh"
+        assert solution["Pch"] > 0, "Chamber pressure should be positive"
+        assert solution["Psub"] > 0, "Vapor pressure should be positive"
+        assert solution["dmdt"] >= 0, "Sublimation rate should be non-negative"
+
+    @pytest.mark.slow
+    def test_solve_with_warmstart(self, standard_vial, standard_product, standard_ht):
+        """Test solving with warmstart initialization."""
+        Lpr0 = functions.Lpr0_FUN(2.0, standard_vial["Ap"], standard_product["cSolid"])
+        Lck = 0.5
+
+        # Create warmstart data with reasonable initial guess
+        warmstart = {
+            "Pch": 0.15,
+            "Tsh": -10.0,
+            "Tsub": -25.0,
+            "Tbot": -20.0,
+            "Psub": 0.5,
+            "dmdt": 0.5,
+            "Kv": 5e-4,
+        }
+
+        model = single_step.create_single_step_model(
+            standard_vial, standard_product, standard_ht, Lpr0, Lck
+        )
+
+        solution = single_step.solve_single_step(
+            model, warmstart_data=warmstart, tee=False
+        )
+
+        # Should solve successfully
+        assert "Pch" in solution
+        assert solution["dmdt"] >= 0
+
+    @pytest.mark.slow
+    def test_optimize_convenience_function(
+        self, standard_vial, standard_product, standard_ht
+    ):
+        """Test the optimize_single_step convenience function."""
+        Lpr0 = functions.Lpr0_FUN(2.0, standard_vial["Ap"], standard_product["cSolid"])
+        Lck = 0.5
+
+        solution = single_step.optimize_single_step(
+            standard_vial, standard_product, standard_ht, Lpr0, Lck, tee=False
+        )
+
+        # Check solution is valid
+        assert "Pch" in solution
+        assert "Tsh" in solution
+        assert solution["dmdt"] >= 0
+
+
+class TestSolutionValidity:
+    """Tests for solution validation utilities."""
+
+    def test_check_valid_solution(self):
+        """Test validation of a physically reasonable solution."""
+        solution = {
+            "Pch": 0.15,
+            "Tsh": -5.0,
+            "Tsub": -25.0,
+            "Tbot": -20.0,
+            "Psub": 0.5,
+            "dmdt": 0.5,
+            "Kv": 5e-4,
+        }
+
+        is_valid, violations = utils.check_solution_validity(solution)
+
+        assert is_valid, f"Valid solution flagged as invalid: {violations}"
+        assert len(violations) == 0
+
+    def test_check_invalid_temperature_ordering(self):
+        """Test detection of invalid temperature ordering."""
+        solution = {
+            "Pch": 0.15,
+            "Tsh": -5.0,
+            "Tsub": -10.0,  # Invalid: Tsub > Tbot
+            "Tbot": -20.0,
+            "Psub": 0.5,
+            "dmdt": 0.5,
+            "Kv": 5e-4,
+        }
+
+        is_valid, violations = utils.check_solution_validity(solution)
+
+        assert not is_valid
+        assert any("Tsub" in v and "Tbot" in v for v in violations)
+
+    def test_check_invalid_driving_force(self):
+        """Test detection of invalid driving force (Psub < Pch)."""
+        solution = {
+            "Pch": 0.5,  # Invalid: Pch > Psub
+            "Tsh": -5.0,
+            "Tsub": -25.0,
+            "Tbot": -20.0,
+            "Psub": 0.3,
+            "dmdt": 0.0,
+            "Kv": 5e-4,
+        }
+
+        is_valid, violations = utils.check_solution_validity(solution)
+
+        assert not is_valid
+        assert any("Psub" in v and "Pch" in v for v in violations)
+
+
+class TestWarmstartUtilities:
+    """Tests for warmstart initialization from scipy."""
+
+    def test_initialize_from_scipy_format(self, standard_vial, standard_product):
+        """Test that warmstart data can be created from scipy output format."""
+        # Create mock scipy output (7 columns)
+        scipy_output = np.array(
+            [
+                [0.0, -25.0, -20.0, -5.0, 150.0, 0.5, 0.3],  # time=0
+                [0.5, -23.0, -18.0, -3.0, 140.0, 0.55, 0.5],  # time=0.5
+            ]
+        )
+
+        Lpr0 = functions.Lpr0_FUN(2.0, standard_vial["Ap"], standard_product["cSolid"])
+
+        warmstart = utils.initialize_from_scipy(
+            scipy_output,
+            time_index=1,
+            vial=standard_vial,
+            product=standard_product,
+            Lpr0=Lpr0,
+        )
+
+        # Check structure
+        assert "Pch" in warmstart
+        assert "Tsh" in warmstart
+        assert "Tsub" in warmstart
+        assert "Tbot" in warmstart
+        assert "Psub" in warmstart
+        assert "dmdt" in warmstart
+        assert "Kv" in warmstart
+
+        # Check values match scipy output (accounting for unit conversions)
+        assert np.isclose(warmstart["Tsub"], -23.0, atol=0.1)
+        assert np.isclose(warmstart["Tbot"], -18.0, atol=0.1)
+        assert np.isclose(warmstart["Tsh"], -3.0, atol=0.1)
+        assert np.isclose(warmstart["Pch"], 0.14, rtol=0.1)  # 140 mTorr → 0.14 Torr
+
+        # Check derived quantities are reasonable
+        assert warmstart["Psub"] > 0
+        assert warmstart["dmdt"] >= 0
+        assert warmstart["Kv"] > 0