fix division by zero error in linear presolve

Author: michaelbynum

pyomo/contrib/solver/ipopt.py

@@ -17,7 +17,11 @@
 
 from pyomo.common import Executable
 from pyomo.common.config import ConfigValue, document_kwargs_from_configdict, ConfigDict
-from pyomo.common.errors import PyomoException, DeveloperError
+from pyomo.common.errors import (
+    PyomoException,
+    DeveloperError,
+    InfeasibleConstraintException,
+)
 from pyomo.common.tempfiles import TempfileManager
 from pyomo.common.timing import HierarchicalTimer
 from pyomo.core.base.var import _GeneralVarData
@@ -314,15 +318,19 @@ def solve(self, model, **kwds):
             ) as row_file, open(basename + '.col', 'w') as col_file:
                 timer.start('write_nl_file')
                 self._writer.config.set_value(config.writer_config)
-                nl_info = self._writer.write(
-                    model,
-                    nl_file,
-                    row_file,
-                    col_file,
-                    symbolic_solver_labels=config.symbolic_solver_labels,
-                )
+                try:
+                    nl_info = self._writer.write(
+                        model,
+                        nl_file,
+                        row_file,
+                        col_file,
+                        symbolic_solver_labels=config.symbolic_solver_labels,
+                    )
+                    proven_infeasible = False
+                except InfeasibleConstraintException:
+                    proven_infeasible = True
                 timer.stop('write_nl_file')
-            if len(nl_info.variables) > 0:
+            if not proven_infeasible and len(nl_info.variables) > 0:
                 # Get a copy of the environment to pass to the subprocess
                 env = os.environ.copy()
                 if nl_info.external_function_libraries:
@@ -361,11 +369,20 @@ def solve(self, model, **kwds):
                     timer.stop('subprocess')
                     # This is the stuff we need to parse to get the iterations
                     # and time
-                    iters, ipopt_time_nofunc, ipopt_time_func, ipopt_total_time = (
-                        self._parse_ipopt_output(ostreams[0])
-                    )
-
-            if len(nl_info.variables) == 0:
+                    (
+                        iters,
+                        ipopt_time_nofunc,
+                        ipopt_time_func,
+                        ipopt_total_time,
+                    ) = self._parse_ipopt_output(ostreams[0])
+
+            if proven_infeasible:
+                results = Results()
+                results.termination_condition = TerminationCondition.provenInfeasible
+                results.solution_loader = SolSolutionLoader(None, None)
+                results.iteration_count = 0
+                results.timing_info.total_seconds = 0
+            elif len(nl_info.variables) == 0:
                 if len(nl_info.eliminated_vars) == 0:
                     results = Results()
                     results.termination_condition = TerminationCondition.emptyModel
@@ -457,7 +474,7 @@ def solve(self, model, **kwds):
                 )
 
         results.solver_configuration = config
-        if len(nl_info.variables) > 0:
+        if not proven_infeasible and len(nl_info.variables) > 0:
             results.solver_log = ostreams[0].getvalue()
 
         # Capture/record end-time / wall-time

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

@@ -1535,6 +1535,40 @@ def test_presolve_with_zero_coef(
         self.assertAlmostEqual(m.y.value, 0)
         self.assertAlmostEqual(m.z.value, 0)
 
+        m.x.setlb(2)
+        res = opt.solve(
+            m, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        self.assertEqual(
+            res.termination_condition, TerminationCondition.provenInfeasible
+        )
+
+        m = pe.ConcreteModel()
+        m.w = pe.Var()
+        m.x = pe.Var()
+        m.y = pe.Var()
+        m.z = pe.Var()
+        m.obj = pe.Objective(expr=m.x**2 + m.y**2 + m.z**2 + m.w**2)
+        m.c1 = pe.Constraint(expr=m.x + m.w == m.y + m.z)
+        m.c2 = pe.Constraint(expr=m.z == -m.y)
+        m.c3 = pe.Constraint(expr=m.x == -m.w)
+
+        res = opt.solve(m)
+        self.assertAlmostEqual(res.incumbent_objective, 0)
+        self.assertAlmostEqual(m.w.value, 0)
+        self.assertAlmostEqual(m.x.value, 0)
+        self.assertAlmostEqual(m.y.value, 0)
+        self.assertAlmostEqual(m.z.value, 0)
+
+        del m.c1
+        m.c1 = pe.Constraint(expr=m.x + m.w == m.y + m.z + 1.5)
+        res = opt.solve(
+            m, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        self.assertEqual(
+            res.termination_condition, TerminationCondition.provenInfeasible
+        )
+
     @parameterized.expand(input=_load_tests(all_solvers))
     def test_scaling(self, name: str, opt_class: Type[SolverBase], use_presolve: bool):
         opt: SolverBase = opt_class()

pyomo/repn/plugins/nl_writer.py

@@ -1781,6 +1781,19 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
                     ):
                         _id, id2 = id2, _id
                         coef, coef2 = coef2, coef
+                if coef == 0 and coef2 == 0:
+                    if abs(expr_info.const) <= TOL:
+                        # constraint is trivially satisfied
+                        # we cannot use this constraint to
+                        # eliminate either variable
+                        eliminated_cons.add(con_id)
+                        continue
+                    else:
+                        # constraint is trivially infeasible
+                        raise InfeasibleConstraintException(
+                            f"Model is infeasible: {expr_info.const} == "
+                            + f"{coef}*{var_map[_id]} + {coef2}*{var_map[id2]}"
+                        )
                 # substituting _id with a*x + b
                 a = -coef2 / coef
                 x = id2
@@ -1799,7 +1812,9 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
                 # Tighten variable bounds
                 lb, ub = var_bounds[_id]
                 if a == 0:
-                    if (lb is not None and lb > b) or (ub is not None and ub < b):
+                    if (lb is not None and lb > b + TOL) or (
+                        ub is not None and ub < b - TOL
+                    ):
                         raise InfeasibleConstraintException(
                             f'Model is infeasible: {lb} <= {var_map[_id]} <= {ub} '
                             + f'and {var_map[_id]} == {a}*{var_map[x]} + {b} cannot '