NLv2: handle presolved independent linear subsystems

pyomo/repn/plugins/nl_writer.py

@@ -1081,9 +1081,37 @@ def write(self, model):
 
         # Update any eliminated variables to point to the (potentially
         # scaled) substituted variables
-        for _id, expr_info in eliminated_vars.items():
+        for _id, expr_info in list(eliminated_vars.items()):
             nl, args, _ = expr_info.compile_repn(visitor)
-            _vmap[_id] = nl.rstrip() % tuple(_vmap[_id] for _id in args)
+            for _i in args:
+                # It is possible that the eliminated variable could
+                # reference another variable that is no longer part of
+                # the model and therefore does not have a _vmap entry.
+                # This can happen when there is an underdetermined
+                # independent linear subsystem and the presolve removed
+                # all the constraints from the subsystem.  Because the
+                # free variables in the subsystem are not referenced
+                # anywhere else in the model, they are not part of the
+                # `variables` list.  Implicitly "fix" it to an arbitrary
+                # valid value from the presolved domain (see #3192).
+                if _i not in _vmap:
+                    lb, ub = var_bounds[_i]
+                    if lb is None:
+                        lb = -inf
+                    if ub is None:
+                        ub = inf
+                    if lb <= 0 <= ub:
+                        val = 0
+                    else:
+                        val = lb if abs(lb) < abs(ub) else ub
+                    eliminated_vars[_i] = AMPLRepn(val, {}, None)
+                    _vmap[_i] = expr_info.compile_repn(visitor)[0]
+                    logger.warning(
+                        "presolve identified an underdetermined independent "
+                        "linear subsystem that was removed from the model.  "
+                        f"Setting '{var_map[_i]}' == {val}"
+                    )
+            _vmap[_id] = nl.rstrip() % tuple(_vmap[_i] for _i in args)
 
         r_lines = [None] * n_cons
         for idx, (con, expr_info, lb, ub) in enumerate(constraints):

pyomo/repn/tests/ampl/test_nlv2.py

@@ -1812,6 +1812,78 @@ def test_presolve_zero_coef(self):
             )
         )
 
+    def test_presolve_independent_subsystem(self):
+        # This is derived from the example in #3192
+        m = ConcreteModel()
+        m.x = Var()
+        m.y = Var()
+        m.z = Var()
+        m.d = Constraint(expr=m.z == m.y)
+        m.c = Constraint(expr=m.y == m.x)
+        m.o = Objective(expr=0)
+
+        ref = """g3 1 1 0       #problem unknown
+ 0 0 1 0 0     #vars, constraints, objectives, ranges, eqns
+ 0 0 0 0 0 0   #nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0   #network constraints: nonlinear, linear
+ 0 0 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 0   #nonzeros in Jacobian, obj. gradient
+ 1 0   #max name lengths: constraints, variables
+ 0 0 0 0 0     #common exprs: b,c,o,c1,o1
+O0 0   #o
+n0
+x0     #initial guess
+r      #0 ranges (rhs's)
+b      #0 bounds (on variables)
+k-1    #intermediate Jacobian column lengths
+"""
+
+        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(),
+            "presolve identified an underdetermined independent linear subsystem "
+            "that was removed from the model.  Setting 'z' == 0\n",
+        )
+
+        self.assertEqual(*nl_diff(ref, OUT.getvalue()))
+
+        m.x.lb = 5.0
+
+        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(),
+            "presolve identified an underdetermined independent linear subsystem "
+            "that was removed from the model.  Setting 'z' == 5.0\n",
+        )
+
+        self.assertEqual(*nl_diff(ref, OUT.getvalue()))
+
+        m.x.lb = -5.0
+        m.z.ub = -2.0
+
+        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(),
+            "presolve identified an underdetermined independent linear subsystem "
+            "that was removed from the model.  Setting 'z' == -2.0\n",
+        )
+
+        self.assertEqual(*nl_diff(ref, OUT.getvalue()))
+
     def test_scaling(self):
         m = pyo.ConcreteModel()
         m.x = pyo.Var(initialize=0)