[WIP] Activating Pyomo linear presolve for IDAES models

idaes/apps/grid_integration/multiperiod/multiperiod.py

@@ -115,7 +115,8 @@ def __init__(
                 self._stochastic_model = False
 
             if solver is None:
-                solver = get_solver()
+                # TODO: fails using linear presolve (non-descript RuntimeError)
+                solver = get_solver(writer_config={"linear_presolve": False})
 
             _logger = logging.getLogger(__name__)
             _logger.setLevel(outlvl)

idaes/config.py

@@ -322,6 +322,106 @@ def _new_idaes_config_block():
         ),
     )
 
+    cfg.declare(
+        "ipopt_v2",
+        pyomo.common.config.ConfigBlock(
+            implicit=False,
+            description="Default config for 'ipopt' solver",
+            doc="Default config for 'ipopt' solver",
+        ),
+    )
+    cfg["ipopt_v2"].declare(
+        "options",
+        pyomo.common.config.ConfigBlock(
+            implicit=True,
+            description="Default solver options for 'ipopt'",
+            doc="Default solver options for 'ipopt' solver",
+        ),
+    )
+
+    cfg["ipopt_v2"]["options"].declare(
+        "nlp_scaling_method",
+        pyomo.common.config.ConfigValue(
+            domain=str,
+            default="gradient-based",
+            description="Ipopt NLP scaling method",
+            doc="Ipopt NLP scaling method",
+        ),
+    )
+
+    cfg["ipopt_v2"]["options"].declare(
+        "tol",
+        pyomo.common.config.ConfigValue(
+            domain=float,
+            default=1e-6,
+            description="Ipopt tol option",
+            doc="Ipopt tol option",
+        ),
+    )
+
+    cfg["ipopt_v2"]["options"].declare(
+        "max_iter",
+        pyomo.common.config.ConfigValue(
+            domain=int,
+            default=200,
+            description="Ipopt max_iter option",
+            doc="Ipopt max_iter option",
+        ),
+    )
+
+    # TODO: Would be nice to make these the defaults, but too many things fail
+    # right now.
+    # TODO: Leaving these as a reminder for the future
+    # cfg["ipopt_v2"]["options"].declare(
+    #     "linear_solver",
+    #     pyomo.common.config.ConfigValue(
+    #         domain=str,
+    #         default="ma57",
+    #         description="Linear solver to be used by IPOPT",
+    #         doc="Linear solver to be used by IPOPT",
+    #     ),
+    # )
+    #
+    # cfg["ipopt_v2"]["options"].declare(
+    #     "ma57_automatic_scaling",
+    #     pyomo.common.config.ConfigValue(
+    #         domain=str,
+    #         default="yes",
+    #         description="Whether to use automatic scaling in MA57",
+    #         doc="Whether to use automatic scaling in MA57. "
+    #         "Valid values: 'yes' (default), 'no'.",
+    #     ),
+    # )
+
+    cfg["ipopt_v2"].declare(
+        "writer_config",
+        pyomo.common.config.ConfigBlock(
+            implicit=True,
+            description="Default writer configuration for 'ipopt'",
+            doc="Default writer configuration for 'ipopt' solver",
+        ),
+    )
+
+    cfg["ipopt_v2"]["writer_config"].declare(
+        "scale_model",
+        pyomo.common.config.ConfigValue(
+            domain=bool,
+            default=False,  # TODO: Change to true once transition complete
+            description="Whether to apply model scaling in writer",
+            doc="Whether to apply model scaling in writer",
+        ),
+    )
+
+    cfg["ipopt_v2"]["writer_config"].declare(
+        "linear_presolve",
+        pyomo.common.config.ConfigValue(
+            domain=bool,
+            default=True,
+            description="Whether to apply linear presolve in writer",
+            doc="Whether to apply linear presolve in writer",
+        ),
+    )
+
     cfg.declare(
         "ipopt_l1",
         pyomo.common.config.ConfigBlock(
@@ -451,7 +551,7 @@ def _new_idaes_config_block():
     cfg.declare(
         "default_solver",
         pyomo.common.config.ConfigValue(
-            default="ipopt",
+            default="ipopt_v2",
             domain=str,
             description="Default solver.  See Pyomo's SolverFactory for details.",
             doc="Default solver.  See Pyomo's SolverFactory for details.",

idaes/core/solvers/config.py

@@ -14,7 +14,10 @@
 # pylint: disable=missing-module-docstring
 # pylint: disable=missing-class-docstring
 
+from copy import deepcopy
+
 from pyomo.environ import SolverFactory
+
 import idaes
 
 
@@ -42,22 +45,31 @@ def __call__(self, *args, **kwargs):
         else:
             name = self.name
             solver = self.solver
+
         if name in idaes.cfg and (
-            idaes.cfg.use_idaes_solver_config
-            or name == "default"
-            or not self.registered
+            idaes.cfg.use_idaes_solver_config or not self.registered
         ):
             for k, v in idaes.cfg[name].items():
                 if k not in kwargs:
-                    kwargs[k] = v
+                    kwargs[k] = deepcopy(v)
                 elif k == "options":
                     # options is in ConfigBlock and in kwargs, treat "options"
                     # special so individual options can have defaults not just
                     # the whole options block
                     for opk, opv in v.items():
                         if opk not in kwargs["options"]:
                             kwargs["options"][opk] = opv
-        return solver(*args, **kwargs)
+                elif k == "writer_config":
+                    # writer_config is in ConfigBlock and in kwargs, treat "writer_config"
+                    # special so individual options can have defaults not just
+                    # the whole options block
+                    for opk, opv in v.items():
+                        if opk not in kwargs["writer_config"]:
+                            kwargs["writer_config"][opk] = opv
+
+        solver_obj = solver(*args, **kwargs)
+
+        return solver_obj
 
 
 def use_idaes_solver_configuration_defaults(b=True):

idaes/core/solvers/get_solver.py

@@ -15,14 +15,16 @@
 This module contains the IDAES get_solver method.
 """
 
+from pyomo.contrib.solver.base import LegacySolverWrapper
+
 import idaes.logger as idaeslog
 import idaes.core.solvers
 
 _log = idaeslog.getLogger(__name__)
 
 
 # Author: Andrew Lee
-def get_solver(solver=None, options=None):
+def get_solver(solver=None, options=None, writer_config=None):
     """
     General method for getting a solver object which defaults to the standard
     IDAES solver (defined in the IDAES configuration).
@@ -32,15 +34,28 @@ def get_solver(solver=None, options=None):
         options: dict of solver options to use, overwrites any settings
                  provided by IDAES configuration. Default = None, use default
                  solver options.
+        writer_config: dict of configuration options for solver writer, overwrites
+                 ny settings provided by IDAES configuration. Default = None, use
+                 default solver options.
 
     Returns:
         A Pyomo solver object
     """
+    # TODO: Deprecate options argument in favour of solver_options to be more consistent
+    # with new Pyomo interfaces?
     if solver is None:
         solver = "default"
     solver_obj = idaes.core.solvers.SolverWrapper(solver, register=False)()
 
-    if options is not None:
-        solver_obj.options.update(options)
+    if isinstance(solver_obj, LegacySolverWrapper):
+        if options is not None:
+            for k, v in options.items():
+                solver_obj.options[k] = v
+        if writer_config is not None:
+            for k, v in writer_config.items():
+                solver_obj.config.writer_config[k] = v
+    else:
+        if options is not None:
+            solver_obj.options.update(options)
 
     return solver_obj

idaes/core/util/initialization.py

@@ -18,6 +18,7 @@
 from pyomo.environ import (
     Block,
     check_optimal_termination,
+    ConcreteModel,
     Constraint,
     value,
 )
@@ -248,7 +249,7 @@ def solve_indexed_blocks(solver, blocks, **kwds):
 
     try:
         # Create a temporary Block
-        tmp = Block(concrete=True)
+        tmp = ConcreteModel(concrete=True)
 
         nBlocks = len(blocks)
 

idaes/core/util/tests/test_model_diagnostics.py

@@ -504,7 +504,8 @@ def model(self):
         m.b.v5.fix(2)
         m.b.v6.fix(0)
 
-        solver = get_solver()
+        # Model is trivially pre-solvable - turn off
+        solver = get_solver(writer_config={"linear_presolve": False})
         solver.solve(m)
 
         return m
@@ -1101,7 +1102,8 @@ def test_collect_numerical_warnings_corrected(self, model):
         m.b.v3.setlb(-5)
         m.b.v5.setub(10)
 
-        solver = get_solver()
+        # Model is trivially pre-solvable - turn off
+        solver = get_solver(writer_config={"linear_presolve": False})
         solver.solve(m)
 
         dt = DiagnosticsToolbox(model=m.b)
@@ -1216,7 +1218,8 @@ def test_assert_no_numerical_warnings(self, model):
         # Fix numerical issues
         m.b.v3.setlb(-5)
 
-        solver = get_solver()
+        # Model is trivially pre-solvable - turn off
+        solver = get_solver(writer_config={"linear_presolve": False})
         solver.solve(m)
 
         dt = DiagnosticsToolbox(model=m.b)

idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py

@@ -758,10 +758,11 @@ def _make_state_vars(self):
                 domain=NonNegativeReals,
                 doc="Total molar flowrate [mol/s]",
                 units=pyunits.mol / pyunits.s,
+                bounds=(1e-8, None),
             )
             self.mole_frac_comp = Var(
                 self.params.component_list,
-                bounds=(0, 1),
+                bounds=(1e-10, 1.1),
                 initialize=1 / len(self.params.component_list),
                 doc="Mixture mole fraction",
             )
@@ -770,12 +771,14 @@ def _make_state_vars(self):
                 domain=NonNegativeReals,
                 doc="State pressure [Pa]",
                 units=pyunits.Pa,
+                bounds=(1e4, None),
             )
             self.temperature = Var(
                 initialize=298.15,
                 domain=NonNegativeReals,
                 doc="State temperature [K]",
                 units=pyunits.K,
+                bounds=(200, None),
             )
         else:
             self.flow_mol_comp = Var(
@@ -784,18 +787,21 @@ def _make_state_vars(self):
                 domain=NonNegativeReals,
                 doc="Component molar flowrate [mol/s]",
                 units=pyunits.mol / pyunits.s,
+                bounds=(1e-8, None),
             )
             self.pressure = Var(
                 initialize=101325,
                 domain=NonNegativeReals,
                 doc="State pressure [Pa]",
                 units=pyunits.Pa,
+                bounds=(1e4, None),
             )
             self.temperature = Var(
                 initialize=298.15,
                 domain=NonNegativeReals,
                 doc="State temperature [K]",
                 units=pyunits.K,
+                bounds=(200, None),
             )
 
     def _make_vars(self):

idaes/models/properties/cubic_eos/cubic_prop_pack.py

@@ -1100,10 +1100,11 @@ def build(self):
             domain=NonNegativeReals,
             doc="Component molar flowrate [mol/s]",
             units=pyunits.mol / pyunits.s,
+            bounds=(1e-8, None),
         )
         self.mole_frac_comp = Var(
             self.params.component_list,
-            bounds=(0, None),
+            bounds=(1e-10, None),
             initialize=1 / len(self.params.component_list),
             doc="Mixture mole fractions [-]",
         )
@@ -1112,12 +1113,14 @@ def build(self):
             domain=NonNegativeReals,
             doc="State pressure [Pa]",
             units=pyunits.Pa,
+            bounds=(1e4, None),
         )
         self.temperature = Var(
             initialize=298.15,
             domain=NonNegativeReals,
             doc="State temperature [K]",
             units=pyunits.K,
+            bounds=(200, None),
         )
 
         # Add supporting variables
@@ -1127,13 +1130,14 @@ def build(self):
             domain=NonNegativeReals,
             doc="Phase molar flow rates [mol/s]",
             units=pyunits.mol / pyunits.s,
+            bounds=(1e-8, None),
         )
 
         self.mole_frac_phase_comp = Var(
             self.params.phase_list,
             self.params.component_list,
             initialize=1 / len(self.params.component_list),
-            bounds=(0, None),
+            bounds=(1e-10, None),
             doc="Phase mole fractions [-]",
         )
 

idaes/models/properties/cubic_eos/tests/test_BT_example.py

@@ -58,6 +58,7 @@ def configure(self):
         self.param_args = {"valid_phase": "Liq"}
         self.prop_args = {}
         self.has_density_terms = True
+        self.skip_initialization_raises_exception_test = True
 
 
 @pytest.mark.skipif(not cubic_roots_available(), reason="Cubic functions not available")
@@ -67,6 +68,7 @@ def configure(self):
         self.param_args = {"valid_phase": "Vap"}
         self.prop_args = {}
         self.has_density_terms = True
+        self.skip_initialization_raises_exception_test = True
 
 
 # -----------------------------------------------------------------------------
@@ -103,7 +105,13 @@ def test_initialization_failure(self):
         m.fs.state[0].temperature.fix(300)
         m.fs.state[0].pressure.fix(1e5)
 
-        m.fs.state[0].p_bound = Constraint(expr=m.fs.state[0].pressure <= 1e4)
+        m.fs.state[0].infeas = Constraint(
+            expr=sum(
+                m.fs.state[0].get_material_flow_terms(p, j) ** 2
+                for p, j in m.fs.state[0].phase_component_set
+            )
+            <= 0
+        )
 
         with pytest.raises(
             InitializationError,