fix division by zero error in linear presolve

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
@@ -72,11 +76,7 @@ def __init__(
             ),
         )
         self.writer_config: ConfigDict = self.declare(
-            'writer_config',
-            ConfigValue(
-                default=NLWriter.CONFIG(),
-                description="Configuration that controls options in the NL writer.",
-            ),
+            'writer_config', NLWriter.CONFIG()
         )
 
 
@@ -314,15 +314,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 +365,17 @@ 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 = (
+                    (iters, ipopt_time_nofunc, ipopt_time_func, ipopt_total_time) = (
                         self._parse_ipopt_output(ostreams[0])
                     )
 
-            if len(nl_info.variables) == 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 +467,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

@@ -1508,6 +1508,71 @@ def test_bug_2(self, name: str, opt_class: Type[SolverBase], use_presolve: bool)
             res = opt.solve(m)
             self.assertAlmostEqual(res.incumbent_objective, -18, 5)
 
+    @parameterized.expand(input=_load_tests(nl_solvers))
+    def test_presolve_with_zero_coef(
+        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 use_presolve:
+            opt.config.writer_config.linear_presolve = True
+        else:
+            opt.config.writer_config.linear_presolve = False
+
+        """
+        when c2 gets presolved out, c1 becomes 
+        x - y + y = 0 which becomes
+        x - 0*y == 0 which is the zero we are testing for
+        """
+        m = pe.ConcreteModel()
+        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.c1 = pe.Constraint(expr=m.x == m.y + m.z + 1.5)
+        m.c2 = pe.Constraint(expr=m.z == -m.y)
+
+        res = opt.solve(m)
+        self.assertAlmostEqual(res.incumbent_objective, 2.25)
+        self.assertAlmostEqual(m.x.value, 1.5)
+        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
+        )
+        if use_presolve:
+            exp = TerminationCondition.provenInfeasible
+        else:
+            exp = TerminationCondition.locallyInfeasible
+        self.assertEqual(res.termination_condition, exp)
+
+        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, exp)
+
     @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

@@ -1760,7 +1760,7 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
                 id2_isdiscrete = var_map[id2].domain.isdiscrete()
                 if var_map[_id].domain.isdiscrete() ^ id2_isdiscrete:
                     # if only one variable is discrete, then we need to
-                    # substiitute out the other
+                    # substitute out the other
                     if id2_isdiscrete:
                         _id, id2 = id2, _id
                         coef, coef2 = coef2, coef
@@ -1820,21 +1820,33 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
                 # appropriately (that expr_info is persisting in the
                 # eliminated_vars dict - and we will use that to
                 # update other linear expressions later.)
+                old_nnz = len(expr_info.linear)
                 c = expr_info.linear.pop(_id, 0)
+                nnz = old_nnz - 1
                 expr_info.const += c * b
                 if x in expr_info.linear:
                     expr_info.linear[x] += c * a
+                    if expr_info.linear[x] == 0:
+                        nnz -= 1
+                        coef = expr_info.linear.pop(x)
                 elif a:
                     expr_info.linear[x] = c * a
                     # replacing _id with x... NNZ is not changing,
                     # but we need to remember that x is now part of
                     # this constraint
                     comp_by_linear_var[x].append((con_id, expr_info))
                     continue
-                # NNZ has been reduced by 1
-                nnz = len(expr_info.linear)
-                _old = lcon_by_linear_nnz[nnz + 1]
+                _old = lcon_by_linear_nnz[old_nnz]
                 if con_id in _old:
+                    if not nnz:
+                        if abs(expr_info.const) > TOL:
+                            # constraint is trivially infeasible
+                            raise InfeasibleConstraintException(
+                                "model contains a trivially infeasible constraint "
+                                f"{expr_info.const} == {coef}*{var_map[x]}"
+                            )
+                        # constraint is trivially feasible
+                        eliminated_cons.add(con_id)
                     lcon_by_linear_nnz[nnz][con_id] = _old.pop(con_id)
             # If variables were replaced by the variable that
             # we are currently eliminating, then we need to update

pyomo/repn/tests/ampl/test_nlv2.py

@@ -1683,6 +1683,131 @@ def test_presolve_named_expressions(self):
 G0 2	#obj
 0 0
 1 0
+""",
+                OUT.getvalue(),
+            )
+        )
+
+    def test_presolve_zero_coef(self):
+        m = ConcreteModel()
+        m.x = Var()
+        m.y = Var()
+        m.z = Var()
+        m.obj = Objective(expr=m.x**2 + m.y**2 + m.z**2)
+        m.c1 = Constraint(expr=m.x == m.y + m.z + 1.5)
+        m.c2 = Constraint(expr=m.z == -m.y)
+
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            nlinfo = nl_writer.NLWriter().write(
+                m, OUT, symbolic_solver_labels=True, linear_presolve=True
+            )
+        self.assertEqual(LOG.getvalue(), "")
+
+        self.assertEqual(nlinfo.eliminated_vars[0], (m.x, 1.5))
+        self.assertIs(nlinfo.eliminated_vars[1][0], m.y)
+        self.assertExpressionsEqual(
+            nlinfo.eliminated_vars[1][1], LinearExpression([-1.0 * m.z])
+        )
+
+        self.assertEqual(
+            *nl_diff(
+                """g3 1 1 0	# problem unknown
+ 1 0 1 0 0	#vars, constraints, objectives, ranges, eqns
+ 0 1 0 0 0 0	#nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0	#network constraints: nonlinear, linear
+ 0 1 0	#nonlinear vars in constraints, objectives, both
+ 0 0 0 1	#linear network variables; functions; arith, flags
+ 0 0 0 0 0	#discrete variables: binary, integer, nonlinear (b,c,o)
+ 0 1	#nonzeros in Jacobian, obj. gradient
+ 3 1	#max name lengths: constraints, variables
+ 0 0 0 0 0	#common exprs: b,c,o,c1,o1
+O0 0	#obj
+o54	#sumlist
+3	#(n)
+o5	#^
+n1.5
+n2
+o5	#^
+o16	#-
+v0	#z
+n2
+o5	#^
+v0	#z
+n2
+x0	#initial guess
+r	#0 ranges (rhs's)
+b	#1 bounds (on variables)
+3	#z
+k0	#intermediate Jacobian column lengths
+G0 1	#obj
+0 0
+""",
+                OUT.getvalue(),
+            )
+        )
+
+        m.c3 = Constraint(expr=m.x == 2)
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            with self.assertRaisesRegex(
+                nl_writer.InfeasibleConstraintException,
+                r"model contains a trivially infeasible constraint 0.5 == 0.0\*y",
+            ):
+                nlinfo = nl_writer.NLWriter().write(
+                    m, OUT, symbolic_solver_labels=True, linear_presolve=True
+                )
+        self.assertEqual(LOG.getvalue(), "")
+
+        m.c1.set_value(m.x >= m.y + m.z + 1.5)
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            nlinfo = nl_writer.NLWriter().write(
+                m, OUT, symbolic_solver_labels=True, linear_presolve=True
+            )
+        self.assertEqual(LOG.getvalue(), "")
+
+        self.assertIs(nlinfo.eliminated_vars[0][0], m.y)
+        self.assertExpressionsEqual(
+            nlinfo.eliminated_vars[0][1], LinearExpression([-1.0 * m.z])
+        )
+        self.assertEqual(nlinfo.eliminated_vars[1], (m.x, 2))
+
+        self.assertEqual(
+            *nl_diff(
+                """g3 1 1 0	# problem unknown
+ 1 1 1 0 0	#vars, constraints, objectives, ranges, eqns
+ 0 1 0 0 0 0	#nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0	#network constraints: nonlinear, linear
+ 0 1 0	#nonlinear vars in constraints, objectives, both
+ 0 0 0 1	#linear network variables; functions; arith, flags
+ 0 0 0 0 0	#discrete variables: binary, integer, nonlinear (b,c,o)
+ 0 1	#nonzeros in Jacobian, obj. gradient
+ 3 1	#max name lengths: constraints, variables
+ 0 0 0 0 0	#common exprs: b,c,o,c1,o1
+C0	#c1
+n0
+O0 0	#obj
+o54	#sumlist
+3	#(n)
+o5	#^
+n2
+n2
+o5	#^
+o16	#-
+v0	#z
+n2
+o5	#^
+v0	#z
+n2
+x0	#initial guess
+r	#1 ranges (rhs's)
+1 0.5	#c1
+b	#1 bounds (on variables)
+3	#z
+k0	#intermediate Jacobian column lengths
+G0 1	#obj
+0 0
 """,
                 OUT.getvalue(),
             )