diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 5194ee0064..4baade5ef0 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -1450,6 +1450,10 @@ def __init__(self, model: _BlockData, **kwargs): self._model, scaled=False, equality_constraints_only=True ) + # Get list of equality constraint and variable names + self._eq_con_list = self.nlp.get_pyomo_equality_constraints() + self._var_list = self.nlp.get_pyomo_variables() + if self.jacobian.shape[0] < 2: raise ValueError( "Model needs at least 2 equality constraints to perform svd_analysis." @@ -1562,10 +1566,6 @@ def display_underdetermined_variables_and_constraints( tol = self.config.size_cutoff_in_singular_vector - # Get list of equality constraint and variable names - eq_con_list = self.nlp.get_pyomo_equality_constraints() - var_list = self.nlp.get_pyomo_variables() - if singular_values is None: singular_values = range(1, len(self.s) + 1) @@ -1588,11 +1588,11 @@ def display_underdetermined_variables_and_constraints( stream.write(f"{TAB}Smallest Singular Value {e}:\n\n") stream.write(f"{2 * TAB}Variables:\n\n") for v in np.where(abs(self.v[:, e - 1]) > tol)[0]: - stream.write(f"{3 * TAB}{var_list[v].name}\n") + stream.write(f"{3 * TAB}{self._var_list[v].name}\n") stream.write(f"\n{2 * TAB}Constraints:\n\n") for c in np.where(abs(self.u[:, e - 1]) > tol)[0]: - stream.write(f"{3 * TAB}{eq_con_list[c].name}\n") + stream.write(f"{3 * TAB}{self._eq_con_list[c].name}\n") stream.write("\n") stream.write("=" * MAX_STR_LENGTH + "\n") @@ -1620,23 +1620,19 @@ def display_constraints_including_variable(self, variable, stream=None): f"object (got {variable})." ) - # Get list of equality constraint and variable names - eq_con_list = self.nlp.get_pyomo_equality_constraints() - var_list = self.nlp.get_pyomo_variables() - # Get index of variable in Jacobian try: - var_idx = var_list.index(variable) - except (ValueError, PyomoException): + var_idx = self.nlp.get_primal_indices([variable])[0] + except (KeyError, PyomoException): raise AttributeError(f"Could not find {variable.name} in model.") - nonzeroes = self.jacobian[:, var_idx].nonzero() + nonzeros = self.jacobian.getcol(var_idx).nonzero() # Build a list of all constraints that include var cons_w_var = [] - for i, r in enumerate(nonzeroes[0]): + for r in nonzeros[0]: cons_w_var.append( - f"{eq_con_list[r].name}: {self.jacobian[(r, nonzeroes[1][i])]}" + f"{self._eq_con_list[r].name}: {self.jacobian[(r, var_idx)]:.3e}" ) # Write the output @@ -1671,23 +1667,20 @@ def display_variables_in_constraint(self, constraint, stream=None): f"object (got {constraint})." ) - # Get list of equality constraint and variable names - eq_con_list = self.nlp.get_pyomo_equality_constraints() - var_list = self.nlp.get_pyomo_variables() - # Get index of variable in Jacobian try: - con_idx = eq_con_list.index(constraint) - except ValueError: + con_idx = self.nlp.get_constraint_indices([constraint])[0] + except KeyError: raise AttributeError(f"Could not find {constraint.name} in model.") - nonzeroes = self.jacobian[con_idx, :].nonzero() + nonzeros = self.jacobian[con_idx, :].nonzero() # Build a list of all vars in constraint vars_in_cons = [] - for i, r in enumerate(nonzeroes[0]): - c = nonzeroes[1][i] - vars_in_cons.append(f"{var_list[c].name}: {self.jacobian[(r, c)]}") + for c in nonzeros[1]: + vars_in_cons.append( + f"{self._var_list[c].name}: {self.jacobian[(con_idx, c)]:.3e}" + ) # Write the output _write_report_section( @@ -2281,8 +2274,8 @@ def __init__(self, block_or_jac, solver=None): self.jac_eq = get_jacobian(self.block, equality_constraints_only=True)[0] # Create a list of equality constraint names - self.eq_con_list = self.nlp.get_pyomo_equality_constraints() - self.var_list = self.nlp.get_pyomo_variables() + self._eq_con_list = self.nlp.get_pyomo_equality_constraints() + self._var_list = self.nlp.get_pyomo_variables() self.candidate_eqns = None @@ -2293,7 +2286,7 @@ def __init__(self, block_or_jac, solver=None): # # TODO: Need to refactor, document, and test support for Jacobian # self.jac_eq = block_or_jac - # self.eq_con_list = None + # self._eq_con_list = None else: raise TypeError("Check the type for 'block_or_jac'") @@ -2813,11 +2806,11 @@ def underdetermined_variables_and_constraints(self, n_calc=1, tol=0.1, dense=Fal ) print("Column: Variable") for i in np.where(abs(self.v[:, n_calc - 1]) > tol)[0]: - print(str(i) + ": " + self.var_list[i].name) + print(str(i) + ": " + self._var_list[i].name) print("") print("Row: Constraint") for i in np.where(abs(self.u[:, n_calc - 1]) > tol)[0]: - print(str(i) + ": " + self.eq_con_list[i].name) + print(str(i) + ": " + self._eq_con_list[i].name) def find_candidate_equations(self, verbose=True, tee=False): """ @@ -2842,7 +2835,7 @@ def find_candidate_equations(self, verbose=True, tee=False): if verbose: print("Solving MILP model...") ce, ds = self._find_candidate_eqs( - self.candidates_milp, self.solver, self.eq_con_list, tee + self.candidates_milp, self.solver, self._eq_con_list, tee ) if ce is not None: @@ -2883,7 +2876,7 @@ def find_irreducible_degenerate_sets(self, verbose=True, tee=False): # Check if equation 'c' is a major element of an IDS ids_ = self._check_candidate_ids( - self.dh_milp, self.solver, c, self.eq_con_list, tee + self.dh_milp, self.solver, c, self._eq_con_list, tee ) if ids_ is not None: diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index de398d5374..3ea6c2e339 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -467,7 +467,7 @@ def model(self): m.b.v2 = Var(units=units.m) m.b.v3 = Var(bounds=(0, 5)) m.b.v4 = Var() - m.b.v5 = Var(bounds=(0, 1)) + m.b.v5 = Var(bounds=(0, 5)) m.b.v6 = Var() m.b.v7 = Var( units=units.m, bounds=(0, 1) @@ -550,7 +550,6 @@ def test_display_variables_at_or_outside_bounds(self, model): The following variable(s) have values at or outside their bounds (tol=0.0E+00): b.v3 (free): value=0.0 bounds=(0, 5) - b.v5 (fixed): value=2 bounds=(0, 1) ==================================================================================== """ @@ -619,7 +618,6 @@ def test_display_variables_near_bounds(self, model): The following variable(s) have values close to their bounds (abs=1.0E-04, rel=1.0E-04): b.v3: value=0.0 bounds=(0, 5) - b.v5: value=2 bounds=(0, 1) b.v7: value=1.0000939326524314e-07 bounds=(0, 1) ==================================================================================== @@ -1009,7 +1007,7 @@ def test_collect_numerical_warnings(self, model): assert len(warnings) == 2 assert "WARNING: 1 Constraint with large residuals (>1.0E-05)" in warnings - assert "WARNING: 2 Variables at or outside bounds (tol=0.0E+00)" in warnings + assert "WARNING: 1 Variable at or outside bounds (tol=0.0E+00)" in warnings assert len(next_steps) == 2 assert "display_constraints_with_large_residuals()" in next_steps @@ -1070,10 +1068,9 @@ def test_collect_numerical_cautions(self, model): dt = DiagnosticsToolbox(model=model.b) cautions = dt._collect_numerical_cautions() - assert len(cautions) == 5 assert ( - "Caution: 3 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)" + "Caution: 2 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)" in cautions ) assert "Caution: 2 Variables with value close to zero (tol=1.0E-08)" in cautions @@ -1133,7 +1130,6 @@ def test_assert_no_numerical_warnings(self, model): # Fix numerical issues m.b.v3.setlb(-5) - m.b.v5.setub(10) solver = get_solver() solver.solve(m) @@ -1198,12 +1194,12 @@ def test_report_numerical_issues(self, model): 2 WARNINGS WARNING: 1 Constraint with large residuals (>1.0E-05) - WARNING: 2 Variables at or outside bounds (tol=0.0E+00) + WARNING: 1 Variable at or outside bounds (tol=0.0E+00) ------------------------------------------------------------------------------------ 5 Cautions - Caution: 3 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04) + Caution: 2 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04) Caution: 2 Variables with value close to zero (tol=1.0E-08) Caution: 1 Variable with extreme value (<1.0E-04 or >1.0E+04) Caution: 1 Variable with None value @@ -1602,8 +1598,8 @@ def test_display_constraints_including_variable(self): expected = """==================================================================================== The following constraints involve v[1]: - c1: 1.0 - c4: 8.0 + c1: 1.000e+00 + c4: 8.000e+00 ==================================================================================== """ @@ -1667,8 +1663,8 @@ def test_display_variables_in_constraint(self): expected = """==================================================================================== The following variables are involved in c1: - v[1]: 1.0 - v[2]: 2.0 + v[1]: 1.000e+00 + v[2]: 2.000e+00 ==================================================================================== """ diff --git a/idaes/core/util/tests/test_utility_minimization.py b/idaes/core/util/tests/test_utility_minimization.py index 696b3b46af..dec7aa8cf3 100644 --- a/idaes/core/util/tests/test_utility_minimization.py +++ b/idaes/core/util/tests/test_utility_minimization.py @@ -127,7 +127,7 @@ def model(self): "state_definition": FTPx, "state_bounds": { "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s), - "temperature": (273.15, 300, 450, pyunits.K), + "temperature": (100, 300, 450, pyunits.K), "pressure": (5e4, 1e5, 1e6, pyunits.Pa), }, "pressure_ref": (1e5, pyunits.Pa), diff --git a/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py b/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py index 672dd413f0..95955efc94 100644 --- a/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py +++ b/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py @@ -80,7 +80,7 @@ def test_T_sweep(self, m): m.fs.obj = Objective(expr=(m.fs.state[1].temperature - 510) ** 2) m.fs.state[1].temperature.setub(600) - for logP in range(8, 13, 1): + for logP in [9.5, 10, 10.5, 11, 11.5, 12]: m.fs.obj.deactivate() m.fs.state[1].flow_mol.fix(100) @@ -115,11 +115,11 @@ def test_P_sweep(self, m): assert check_optimal_termination(results) while m.fs.state[1].pressure.value <= 1e6: - m.fs.state[1].pressure.value = m.fs.state[1].pressure.value + 1e5 results = solver.solve(m) assert check_optimal_termination(results) - print(T, m.fs.state[1].pressure.value) + + m.fs.state[1].pressure.value = m.fs.state[1].pressure.value + 1e5 @pytest.mark.component def test_T350_P1_x5(self, m): diff --git a/idaes/models/unit_models/tests/test_equilibrium_reactor.py b/idaes/models/unit_models/tests/test_equilibrium_reactor.py index 495f70b6ae..dd1b07047d 100644 --- a/idaes/models/unit_models/tests/test_equilibrium_reactor.py +++ b/idaes/models/unit_models/tests/test_equilibrium_reactor.py @@ -253,13 +253,6 @@ def test_get_performance_contents(self, sapon): } } - @pytest.mark.component - def test_initialization_error(self, sapon): - sapon.fs.unit.outlet.pressure[0].fix(1) - - with pytest.raises(InitializationError): - sapon.fs.unit.initialize() - class TestInitializers: @pytest.fixture diff --git a/idaes/models/unit_models/tests/test_hx_ntu.py b/idaes/models/unit_models/tests/test_hx_ntu.py index 8bf5ec0f7b..e382246652 100644 --- a/idaes/models/unit_models/tests/test_hx_ntu.py +++ b/idaes/models/unit_models/tests/test_hx_ntu.py @@ -505,16 +505,6 @@ def test_conservation(self, model): <= 1e-6 ) - @pytest.mark.component - def test_initialization_error(self, model): - model.fs.unit.hot_side_outlet.pressure[0].fix(1) - - with pytest.raises(InitializationError): - model.fs.unit.initialize() - - # Revert DoF change to avoid contaminating subsequent tests - model.fs.unit.hot_side_outlet.pressure[0].unfix() - class TestInitializers(object): @pytest.fixture(scope="class") diff --git a/idaes/models_extra/column_models/plate_heat_exchanger.py b/idaes/models_extra/column_models/plate_heat_exchanger.py index 0fc4dc98b0..15a52f6381 100644 --- a/idaes/models_extra/column_models/plate_heat_exchanger.py +++ b/idaes/models_extra/column_models/plate_heat_exchanger.py @@ -46,6 +46,7 @@ ) from pyomo.common.config import ConfigValue, In, Integer from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.common.deprecation import deprecated # Import IDAES cores from idaes.core import declare_process_block_class @@ -61,6 +62,13 @@ _log = idaeslog.getLogger(__name__) +@deprecated( + "The Plate Heat Exchanger (PHE) model is known to be affected by " + "issues causing it to fail to solve on certain platforms starting with Pyomo v6.7.0. " + "This might cause the model to be removed in an upcoming IDAES release if these failures are not resolved. " + "For more information, see IDAES/idaes-pse#1294", + version="2.3.0", +) @declare_process_block_class("PlateHeatExchanger") class PlateHeatExchangerData(HeatExchangerNTUData): """Plate Heat Exchanger(PHE) Unit Model.""" diff --git a/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py b/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py index 1d05c9e0b7..ce3983c753 100644 --- a/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py +++ b/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py @@ -65,6 +65,12 @@ def test_config(): assert len(m.fs.unit.config) == 9 +workaround_for_1294 = pytest.mark.xfail( + reason="These tests fail with Pyomo 6.7.0. See IDAES/idaes-pse#1294 for details", + strict=False, # the failures only occur for certain platforms, e.g. Windows on GHA +) + + # ----------------------------------------------------------------------------- class TestPHE(object): @pytest.fixture(scope="class") @@ -166,6 +172,7 @@ def test_initialize(self, phe): phe, duty=(245000, pyunits.W), optarg={"bound_push": 1e-8, "mu_init": 1e-8} ) + @workaround_for_1294 @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component @@ -175,6 +182,7 @@ def test_solve(self, phe): # Check for optimal solution assert check_optimal_termination(results) + @workaround_for_1294 @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component @@ -214,6 +222,7 @@ def test_solution(self, phe): phe.fs.unit.cold_side_outlet.temperature[0] ) + @workaround_for_1294 @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @pytest.mark.component diff --git a/idaes/models_extra/temperature_swing_adsorption/fixed_bed_tsa0d.py b/idaes/models_extra/temperature_swing_adsorption/fixed_bed_tsa0d.py index 62d14c6aad..62f0a22a8b 100644 --- a/idaes/models_extra/temperature_swing_adsorption/fixed_bed_tsa0d.py +++ b/idaes/models_extra/temperature_swing_adsorption/fixed_bed_tsa0d.py @@ -1167,29 +1167,35 @@ def pressure_n2_rich_stream_eq(b, t): self.flow_mol_h2o_o2_stream = Var( self.flowsheet().time, self.component_list - self.isotherm_components, + initialize=1, units=units.mol / units.s, doc="Component mole flow rate at H2O+O2 stream", ) self.temperature_h2o_o2_stream = Var( - self.flowsheet().time, units=units.K, doc="Temperature at H2O+O2 stream" + self.flowsheet().time, + initialize=300, + units=units.K, + doc="Temperature at H2O+O2 stream", ) self.pressure_h2o_o2_stream = Var( - self.flowsheet().time, units=units.Pa, doc="Pressure at H2O+O2 stream" + self.flowsheet().time, + initialize=101325, + units=units.Pa, + doc="Pressure at H2O+O2 stream", ) - # create empty port - p = Port(noruleinit=True, doc="Outlet port for H2O+O2 stream") - # add port object as an attribute to model - setattr(self, "h2o_o2_stream", p) # dictionary containing state of outlet h2o_o2_stream h2o_o2_stream_dict = { "flow_mol_comp": self.flow_mol_h2o_o2_stream, "temperature": self.temperature_h2o_o2_stream, "pressure": self.pressure_h2o_o2_stream, } - # populate port and map names to actual variables as defined - for key, item in h2o_o2_stream_dict.items(): - p.add(item, name=key) + # create port + self.h2o_o2_stream = Port( + noruleinit=True, + initialize=h2o_o2_stream_dict, + doc="Outlet port for H2O+O2 stream", + ) # add constraints to populate h2o_o2_stream @self.Constraint( @@ -1632,7 +1638,9 @@ def _add_pressurization_step(self): units=units.dimensionless, doc="Mole fraction in pressurization step", ) - self.pressurization.time = Var(units=units.s, doc="Time of pressurization step") + self.pressurization.time = Var( + initialize=1e2, units=units.s, doc="Time of pressurization step" + ) self.pressurization.loading = Var( self.isotherm_components, units=units.mol / units.kg, @@ -1758,7 +1766,9 @@ def _add_adsorption_step(self): self.adsorption = SkeletonUnitModel() # variables - self.adsorption.time = Var(units=units.s, doc="Time of adsorption step") + self.adsorption.time = Var( + initialize=1e2, units=units.s, doc="Time of adsorption step" + ) self.adsorption.loading = Var( self.isotherm_components, units=units.mol / units.kg, @@ -2843,7 +2853,7 @@ def _get_stream_table_contents(self, time_point=0): "Inlet": "inlet", "CO2 Rich Stream": "co2_rich_stream", "N2 Rich Stream": "n2_rich_stream", - "H20 Stream": "h2o_o2_stream", + "H2O Stream": "h2o_o2_stream", }.items(): port_obj = getattr(self, v) diff --git a/idaes/models_extra/temperature_swing_adsorption/tests/test_fixed_bed_tsa0d.py b/idaes/models_extra/temperature_swing_adsorption/tests/test_fixed_bed_tsa0d.py index e0a95d2e8d..93ab2ed8cd 100644 --- a/idaes/models_extra/temperature_swing_adsorption/tests/test_fixed_bed_tsa0d.py +++ b/idaes/models_extra/temperature_swing_adsorption/tests/test_fixed_bed_tsa0d.py @@ -17,6 +17,7 @@ __author__ = "Alex Noring" import pytest +import numpy as np from pyomo.environ import ( check_optimal_termination, @@ -161,11 +162,28 @@ def test_solution(self, model): # air velocity assert pytest.approx(2.44178, abs=1e-3) == value(model.fs.unit.velocity_in) - @pytest.mark.ui @pytest.mark.unit def test_report(self, model): - model.fs.unit.report() - tsa_summary(model.fs.unit) + stream_table_df = model.fs.unit._get_stream_table_contents() + expected_columns = ["Inlet", "CO2 Rich Stream", "N2 Rich Stream", "H2O Stream"] + for column in expected_columns: + assert column in stream_table_df.columns + assert not np.isnan(stream_table_df[column]["temperature"]) + assert not np.isnan(stream_table_df[column]["pressure"]) + + performance_dict = model.fs.unit._get_performance_contents() + for k in performance_dict["vars"]: + assert performance_dict["vars"][k] is not None + + @pytest.mark.unit + def test_summary(self, model): + summary_df = tsa_summary(model.fs.unit) + assert "Adsorption temperature [K]" in summary_df.index + assert "Cycle time [h]" in summary_df.index + assert "Pressure drop [Pa]" in summary_df.index + + @pytest.mark.unit + def test_plotting(self, model): plot_tsa_profiles(model.fs.unit) diff --git a/idaes/models_extra/temperature_swing_adsorption/util.py b/idaes/models_extra/temperature_swing_adsorption/util.py index 6e84988a3d..8e3d4d3b0e 100644 --- a/idaes/models_extra/temperature_swing_adsorption/util.py +++ b/idaes/models_extra/temperature_swing_adsorption/util.py @@ -63,6 +63,8 @@ def tsa_summary(tsa, stream=stdout, export=False): stream.write(textwrap.indent(stream_table_dataframe_to_string(df), " " * 4)) stream.write("\n" + "=" * 84 + "\n") + return df + def plot_tsa_profiles(tsa): """ diff --git a/setup.py b/setup.py index 6887335c28..558f7377b6 100644 --- a/setup.py +++ b/setup.py @@ -127,7 +127,7 @@ def __getitem__(self, key): # Put abstract (non-versioned) deps here. # Concrete dependencies go in requirements[-dev].txt install_requires=[ - "pyomo @ git+https://github.com/IDAES/Pyomo@6.7.0.idaes.2023.11.6", + "pyomo >= 6.7.0", "pint", # required to use Pyomo units "networkx", # required to use Pyomo network "numpy",