Merge pull request #3673 from jsiirola/nl-scaling-named-expr

NL writer: resolve model scaling bug writing `Expression` objects

Author: blnicho

pyomo/repn/plugins/nl_writer.py

@@ -609,6 +609,7 @@ def write(self, model):
             del suffix_data['scaling_factor']
         else:
             scaling_factor = _NoScalingFactor()
+            scaling_cache = None
         scale_model = scaling_factor.scale
 
         timer.toc("Collected suffixes", level=logging.DEBUG)
@@ -1360,7 +1361,7 @@ def write(self, model):
                 # Note: checking target_expr == 0 is equivalent to
                 # testing "(_con_id is not None and _obj_id is not None)
                 # or _con_id == 0 or _obj_id == 0"
-                self._write_v_line(_id, 0)
+                self._write_v_line(_id, 0, scale_model, scaling_cache)
             else:
                 if target_expr not in single_use_subexpressions:
                     single_use_subexpressions[target_expr] = []
@@ -1395,7 +1396,7 @@ def write(self, model):
                 break
             if single_use_subexpressions:
                 for _id in single_use_subexpressions.get(id(info[0]), ()):
-                    self._write_v_line(_id, row_idx + 1)
+                    self._write_v_line(_id, row_idx + 1, scale_model, scaling_cache)
             ostream.write(f'C{row_idx}{row_comments[row_idx]}\n')
             self._write_nl_expression(info[1], False)
 
@@ -1408,7 +1409,9 @@ def write(self, model):
                     # Note that "Writing .nl files" (2005) is incorrectly
                     # missing the "+ 1" in the description of V lines
                     # appearing in only Objectives (bottom of page 9).
-                    self._write_v_line(_id, n_cons + n_lcons + obj_idx + 1)
+                    self._write_v_line(
+                        _id, n_cons + n_lcons + obj_idx + 1, scale_model, scaling_cache
+                    )
             lbl = row_comments[n_cons + obj_idx]
             sense = 0 if info[0].sense == minimize else 1
             ostream.write(f'O{obj_idx} {sense}{lbl}\n')
@@ -1982,7 +1985,7 @@ def _write_nl_expression(self, repn, include_const):
         else:
             self.ostream.write(self.template.const % 0)
 
-    def _write_v_line(self, expr_id, k):
+    def _write_v_line(self, expr_id, k, scale_model, scaling_cache):
         ostream = self.ostream
         column_order = self.column_order
         info = self.subexpression_cache[expr_id]
@@ -1994,10 +1997,14 @@ def _write_v_line(self, expr_id, k):
         # Do NOT write out 0 coefficients here: doing so fouls up the
         # ASL's logic for calculating derivatives, leading to 'nan' in
         # the Hessian results.
-        linear = dict(item for item in info[1].linear.items() if item[1])
+        linear = info[1].linear
+        linear_ids = list(_id for _id, coef in linear.items() if coef)
+        if scale_model:
+            for _id in linear_ids:
+                linear[_id] /= scaling_cache[_id]
         #
-        ostream.write(f'V{self.next_V_line_id} {len(linear)} {k}{lbl}\n')
-        for _id in sorted(linear, key=column_order.__getitem__):
+        ostream.write(f'V{self.next_V_line_id} {len(linear_ids)} {k}{lbl}\n')
+        for _id in sorted(linear_ids, key=column_order.__getitem__):
             ostream.write(f'{column_order[_id]} {linear[_id]!s}\n')
         self._write_nl_expression(info[1], True)
         self.next_V_line_id += 1

pyomo/repn/tests/ampl/test_nlv2.py

@@ -2117,6 +2117,92 @@ def test_scaling(self):
         # Debugging: this diffs the unscaled & scaled models
         # self.assertEqual(*nl_diff(nl1, nl2))
 
+    def test_scaling_named_expressions(self):
+        m = pyo.ConcreteModel()
+
+        m.x = pyo.Var(initialize=1 / 2)
+        m.y = pyo.Var(initialize=1 / 4)
+
+        m.lin_expr = pyo.Expression(expr=m.x + m.y + m.x**3 + m.y**5)
+        m.eqn1 = pyo.Constraint(expr=m.lin_expr**2 == 2)
+        m.eqn2 = pyo.Constraint(expr=m.x**2 + m.y == 3)
+
+        m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT)
+        m.scaling_factor[m.x] = 2
+        m.scaling_factor[m.y] = 4
+        m.scaling_factor[m.eqn1] = 0.5
+
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            nlinfo = nl_writer.NLWriter().write(
+                m, OUT, scale_model=True, linear_presolve=False
+            )
+        self.assertEqual(LOG.getvalue(), "")
+
+        nl2 = OUT.getvalue()
+
+        self.assertEqual(
+            *nl_diff(
+                """g3 1 1 0       #problem unknown
+ 2 2 0 0 2     #vars, constraints, objectives, ranges, eqns
+ 2 0 0 0 0 0   #nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0   #network constraints: nonlinear, linear
+ 2 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)
+ 4 0   #nonzeros in Jacobian, obj. gradient
+ 0 0   #max name lengths: constraints, variables
+ 0 0 0 2 0     #common exprs: b,c,o,c1,o1
+V2 0 1
+o0
+o5
+o3
+v0
+n2
+n3
+o5
+o3
+v1
+n4
+n5
+V3 2 1
+0 0.5
+1 0.25
+v2
+C0
+o2
+n0.5
+o5
+v3
+n2
+C1
+o5
+o3
+v0
+n2
+n2
+x2
+0 1
+1 1
+r
+4 1
+4 3
+b
+3
+3
+k1
+2
+J0 2
+0 0
+1 0
+J1 2
+0 0
+1 0.25
+""",
+                nl2,
+            )
+        )
+
     def test_named_expressions(self):
         # This tests an error possibly reported by #2810
         m = ConcreteModel()