Add tests for trivial constraints in pyomo/contrib/solver

pyomo/contrib/observer/model_observer.py

@@ -614,7 +614,7 @@ def _check_for_var_changes(self):
                 vars_to_update.append(v)
             elif _domain_interval != v.domain.get_interval():
                 vars_to_update.append(v)
-            elif v.value != _value:
+            elif v.fixed and v.value != _value:
                 vars_to_update.append(v)
         cons_to_update = list(cons_to_update.keys())
         return vars_to_update, cons_to_update, update_obj

pyomo/contrib/solver/common/solution_loader.py

@@ -17,6 +17,7 @@
 from pyomo.core.base.var import VarData
 from pyomo.core.staleflag import StaleFlagManager
 from pyomo.core.base.suffix import Suffix
+from .util import NoSolutionError
 
 
 def load_import_suffixes(
@@ -187,6 +188,43 @@ def load_import_suffixes(self, solution_id=None):
         return NotImplemented
 
 
+class NoSolutionSolutionLoader(SolutionLoaderBase):
+    def __init__(self) -> None:
+        pass
+
+    def get_solution_ids(self) -> List[Any]:
+        return []
+
+    def get_number_of_solutions(self) -> int:
+        return 0
+
+    def load_solution(self, solution_id=None):
+        raise NoSolutionError()
+
+    def load_vars(
+        self, vars_to_load: Sequence[VarData] | None = None, solution_id=None
+    ) -> None:
+        raise NoSolutionError()
+
+    def get_vars(
+        self, vars_to_load: Sequence[VarData] | None = None, solution_id=None
+    ) -> Mapping[VarData, float]:
+        raise NoSolutionError()
+
+    def get_duals(
+        self, cons_to_load: Sequence[ConstraintData] | None = None, solution_id=None
+    ) -> Dict[ConstraintData, float]:
+        raise NoSolutionError()
+
+    def get_reduced_costs(
+        self, vars_to_load: Sequence[VarData] | None = None, solution_id=None
+    ) -> Mapping[VarData, float]:
+        raise NoSolutionError()
+
+    def load_import_suffixes(self, solution_id=None):
+        raise NoSolutionError()
+
+
 class PersistentSolutionLoader(SolutionLoaderBase):
     """
     Loader for persistent solvers

pyomo/contrib/solver/solvers/gurobi/gurobi_direct_base.py

@@ -18,7 +18,7 @@
 from pyomo.common.config import ConfigValue
 from pyomo.common.dependencies import attempt_import
 from pyomo.common.enums import ObjectiveSense
-from pyomo.common.errors import ApplicationError
+from pyomo.common.errors import ApplicationError, InfeasibleConstraintException
 from pyomo.common.shutdown import python_is_shutting_down
 from pyomo.common.tee import capture_output, TeeStream
 from pyomo.common.timing import HierarchicalTimer
@@ -33,6 +33,7 @@
     NoReducedCostsError,
     NoSolutionError,
 )
+from pyomo.contrib.solver.common.solution_loader import NoSolutionSolutionLoader
 from pyomo.contrib.solver.common.results import (
     Results,
     SolutionStatus,
@@ -353,6 +354,8 @@ def solve(self, model, **kwds) -> Results:
             res = self._postsolve(
                 grb_model=gurobi_model, solution_loader=solution_loader, has_obj=has_obj
             )
+        except InfeasibleConstraintException:
+            res = self._get_infeasible_results()
         finally:
             os.chdir(orig_cwd)
 
@@ -390,6 +393,24 @@ def _get_tc_map(self):
             }
         return GurobiDirectBase._tc_map
 
+    def _get_infeasible_results(self):
+        res = Results()
+        res.solution_loader = NoSolutionSolutionLoader()
+        res.solution_status = SolutionStatus.noSolution
+        res.termination_condition = TerminationCondition.provenInfeasible
+        res.incumbent_objective = None
+        res.objective_bound = None
+        res.iteration_count = None
+        res.timing_info.gurobi_time = None
+        res.solver_config = self.config
+        res.solver_name = self.name
+        res.solver_version = self.version()
+        if self.config.raise_exception_on_nonoptimal_result:
+            raise NoOptimalSolutionError()
+        if self.config.load_solutions:
+            raise NoFeasibleSolutionError()
+        return res
+
     def _postsolve(self, grb_model, solution_loader, has_obj):
         status = grb_model.Status
 

pyomo/contrib/solver/solvers/gurobi/gurobi_persistent.py

@@ -17,6 +17,7 @@
 from pyomo.common.collections import ComponentSet, OrderedSet
 from pyomo.common.shutdown import python_is_shutting_down
 from pyomo.common.timing import HierarchicalTimer
+from pyomo.common.errors import InfeasibleConstraintException
 from pyomo.core.base.objective import ObjectiveData
 from pyomo.core.kernel.objective import minimize, maximize
 from pyomo.core.base.var import VarData
@@ -487,7 +488,9 @@ def _get_expr_from_pyomo_repn(self, repn):
             vlist = [self._pyomo_var_to_solver_var_map[id(v)] for v in repn.linear_vars]
             new_expr = gurobipy.LinExpr(coef_list, vlist)
         else:
-            new_expr = 0.0
+            # this can't just be zero in case the constraint is a
+            # trivial one
+            new_expr = gurobipy.LinExpr()
 
         if len(repn.quadratic_vars) > 0:
             missing_vars = {}
@@ -737,6 +740,7 @@ def _add_constraints(self, cons: List[ConstraintData]):
         for ndx, con in enumerate(cons):
             lb, body, ub = con.to_bounded_expression(evaluate_bounds=False)
             repn = generate_standard_repn(body, quadratic=True, compute_values=False)
+
             if len(repn.quadratic_vars) > 0:
                 self._quadratic_cons.add(con)
             else:

pyomo/contrib/solver/tests/solvers/test_solvers.py

@@ -35,6 +35,11 @@
     GurobiDirectQuadratic,
     GurobiPersistent,
 )
+from pyomo.contrib.solver.common.util import (
+    NoSolutionError,
+    NoFeasibleSolutionError,
+    NoOptimalSolutionError,
+)
 from pyomo.contrib.solver.solvers.highs import Highs
 from pyomo.core.expr.numeric_expr import LinearExpression
 from pyomo.core.expr.compare import assertExpressionsEqual
@@ -1053,6 +1058,58 @@ def test_results_infeasible(
             ):
                 res.solution_loader.get_reduced_costs()
 
+    @parameterized.expand(input=_load_tests(all_solvers))
+    def test_trivial_constraints(
+        self, name: str, opt_class: Type[SolverBase], use_presolve: bool
+    ):
+        opt: SolverBase = opt_class()
+        if not opt.available():
+            raise unittest.SkipTest(f'Solver {opt.name} not available.')
+        if any(name.startswith(i) for i in nl_solvers_set):
+            if use_presolve:
+                opt.config.writer_config.linear_presolve = True
+            else:
+                opt.config.writer_config.linear_presolve = False
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var()
+        m.y = pyo.Var()
+        m.obj = pyo.Objective(expr=m.y)
+        m.c1 = pyo.Constraint(expr=m.y >= m.x)
+        m.c2 = pyo.Constraint(expr=m.y >= -m.x)
+        m.c3 = pyo.Constraint(expr=m.x >= 0)
+
+        res = opt.solve(m)
+        self.assertAlmostEqual(m.x.value, 0)
+        self.assertAlmostEqual(m.y.value, 0)
+
+        m.x.fix(1)
+        opt.config.tee = True
+        res = opt.solve(m)
+        self.assertAlmostEqual(m.x.value, 1)
+        self.assertAlmostEqual(m.y.value, 1)
+
+        m.x.fix(-1)
+        with self.assertRaises(NoOptimalSolutionError):
+            res = opt.solve(m)
+
+        opt.config.load_solutions = False
+        opt.config.raise_exception_on_nonoptimal_result = False
+        res = opt.solve(m)
+        self.assertNotEqual(res.solution_status, SolutionStatus.optimal)
+        if isinstance(opt, Ipopt):
+            acceptable_termination_conditions = {
+                TerminationCondition.locallyInfeasible,
+                TerminationCondition.unbounded,
+                TerminationCondition.provenInfeasible,
+            }
+        else:
+            acceptable_termination_conditions = {
+                TerminationCondition.provenInfeasible,
+                TerminationCondition.infeasibleOrUnbounded,
+            }
+        self.assertIn(res.termination_condition, acceptable_termination_conditions)
+        self.assertIsNone(res.incumbent_objective)
+
     @parameterized.expand(input=_load_tests(all_solvers))
     def test_duals(self, name: str, opt_class: Type[SolverBase], use_presolve: bool):
         opt: SolverBase = opt_class()

pyomo/repn/plugins/standard_form.py

@@ -20,6 +20,7 @@
     InEnum,
     document_kwargs_from_configdict,
 )
+from pyomo.common.errors import InfeasibleConstraintException
 from pyomo.common.dependencies import scipy, numpy as np
 from pyomo.common.enums import ObjectiveSense
 from pyomo.common.gc_manager import PauseGC
@@ -462,7 +463,7 @@ def write(self, model):
                 # TODO: add a (configurable) feasibility tolerance
                 if (lb is None or lb <= offset) and (ub is None or ub >= offset):
                     continue
-                raise InfeasibleError(
+                raise InfeasibleConstraintException(
                     f"model contains a trivially infeasible constraint, '{con.name}'"
                 )