fix division by zero error in linear presolve

[michaelbynum]

Summary/Motivation:

In [1]: import pyomo.environ as pe
   ...: from pyomo.repn.plugins.nl_writer import NLWriter
   ...: import io
   ...: 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)
   ...: writer = NLWriter()
   ...: out = io.StringIO()
   ...: writer.write(m, out)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 13
     11 writer = NLWriter()
     12 out = io.StringIO()
---> 13 writer.write(m, out)

File /.../pyomo/repn/plugins/nl_writer.py:409, in NLWriter.write(self, model, ostream, rowstream, colstream, **options)
    405 # Pause the GC, as the walker that generates the compiled NL
    406 # representation generates (and disposes of) a large number of
    407 # small objects.
    408 with _NLWriter_impl(ostream, rowstream, colstream, config) as impl:
--> 409     return impl.write(model)

File /.../pyomo/repn/plugins/nl_writer.py:775, in _NLWriter_impl.write(self, model)
    769 # This may fetch more bounds than needed, but only in the cases
    770 # where variables were completely eliminated while walking the
    771 # expressions, or when users provide superfluous variables in
    772 # the column ordering.
    773 var_bounds = {_id: v.bounds for _id, v in var_map.items()}
--> 775 eliminated_cons, eliminated_vars = self._linear_presolve(
    776     comp_by_linear_var, lcon_by_linear_nnz, var_bounds
    777 )
    778 del comp_by_linear_var
    779 del lcon_by_linear_nnz

File /.../pyomo/repn/plugins/nl_writer.py:1771, in _NLWriter_impl._linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds)
   1766         coef, coef2 = coef2, coef
   1767 else:
   1768     # In an attempt to improve numerical stability, we will
   1769     # solve for (and substitute out) the variable with the
   1770     # coefficient closer to +/-1)
-> 1771     log_coef = _log10(abs(coef))
   1772     log_coef2 = _log10(abs(coef2))
   1773     if abs(log_coef2) < abs(log_coef) or (
   1774         log_coef2 == -log_coef and log_coef2 < log_coef
   1775     ):

ValueError: math domain error

Changes proposed in this PR:

• Check for division by zero in NL writer's linear presolve

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

1. I agree my contributions are submitted under the BSD license.
2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

[jsiirola]

Will you still get a divide-by-0 error if both coefficients are 0? Maybe something like:

m.c1 = pe.Constraint(expr=m.x + m.w == m.y + m.z + 1.5)
m.c2 = pe.Constraint(expr=m.z == -m.y)
m.c3 = pe.Constraint(expr=m.w == -m.x)

Maybe the right answer is to filter out 0's when we are performing the substitution?

[michaelbynum]

If both coeffiicents are zero, the constraint should simplify to trivially satisfied (if b is 0) or trivially infeasible (if b is not 0). Let me test for that real quick.

[michaelbynum]

Okay, take a look at this version.

[jsiirola]

I think this can work, but I wonder if there is a simpler implementation: (I think) the only way to get 0 coefficients is through previous substitutions. If that is the case, then, you could revert all the edits to this file and only change (starting around @ L1853):

                 # 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)
+                        if not nnz:
+                            if abs(expr_info.const) > TOL:
+                                # constraint is trivially infeasible
+                                raise InfeasibleConstraintException(
+                                    "model contains a trivially infeasible constrint "
+                                    f"{expr_info.const} == {coef}*{var_map[x]}"
+                                )
+                            # constraint is trivially feasible
+                            eliminated_cons.add(con_id)
                 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:
                     lcon_by_linear_nnz[nnz][con_id] = _old.pop(con_id)

[michaelbynum]

That worked! I'll push shortly.

[michaelbynum]

Okay, third time is the charm! Check out this version.

P.s. git is the best.

[mrmundt]

I'm not part of this conversation, but I had to +1 this. Git is the best.

[jsiirola]

Ick. I think I lied to you. The block beginning if nnz == 0: should be moved farther down:

@@ -1829,15 +1829,6 @@ class _NLWriter_impl(object):
                     if expr_info.linear[x] == 0:
                         nnz -= 1
                         coef = expr_info.linear.pop(x)
-                        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)
                 elif a:
                     expr_info.linear[x] = c * a
                     # replacing _id with x... NNZ is not changing,
@@ -1847,6 +1838,15 @@ class _NLWriter_impl(object):
                     continue
                 _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 constrint "
+                                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

The issue is in the original location, it will check the linear NNZ for all constraints (including inequalities and nonlinear constraints). In the new location, it will only check for the linear equality constraints (which is what we want).

[michaelbynum]

Done.