Preserve affine GDPopt expressions with mutable-parameter coefficients

pyomo/contrib/gdpopt/branch_and_bound.py

@@ -34,16 +34,15 @@
 from pyomo.contrib.gdpopt.nlp_initialization import restore_vars_to_original_values
 from pyomo.contrib.gdpopt.util import (
     copy_var_list_values,
-    SuppressInfeasibleWarning,
     get_main_elapsed_time,
+    is_affine,
+    SuppressInfeasibleWarning,
 )
 from pyomo.contrib.satsolver.satsolver import satisfiable
 from pyomo.core import minimize, Suffix, Constraint, TransformationFactory
 from pyomo.opt import SolverFactory, SolverStatus
 from pyomo.opt import TerminationCondition as tc
 
-_linear_degrees = {1, 0}
-
 # Data tuple for each node that also functions as the sort key.
 # Therefore, ordering of the arguments below matters.
 BBNodeData = namedtuple(
@@ -166,7 +165,7 @@ def _solve_gdp(self, model, config):
                     for constr in disjunct.component_data_objects(
                         Constraint, active=True
                     )
-                    if constr.body.polynomial_degree() not in _linear_degrees
+                    if not is_affine(constr.body)
                 ]
                 for constraint in nonlinear_constraints_in_disjunct:
                     constraint.deactivate()

pyomo/contrib/gdpopt/create_oa_subproblems.py

@@ -21,6 +21,7 @@
 from pyomo.contrib.gdpopt.discrete_problem_initialize import valid_init_strategies
 from pyomo.contrib.gdpopt.util import (
     get_main_elapsed_time,
+    is_affine,
     move_nonlinear_objective_to_constraints,
 )
 from pyomo.gdp.disjunct import Disjunct, Disjunction
@@ -75,7 +76,7 @@ def initialize_discrete_problem(util_block, subprob_util_block, config, solver):
     for c in discrete.component_data_objects(
         Constraint, active=True, descend_into=(Block, Disjunct)
     ):
-        if c.body.polynomial_degree() not in (1, 0):
+        if not is_affine(c.body):
             c.deactivate()
 
     # Transform to a MILP

pyomo/contrib/gdpopt/discrete_problem_initialize.py

@@ -18,7 +18,7 @@
 
 from pyomo.contrib.gdpopt.cut_generation import add_no_good_cut
 from pyomo.contrib.gdpopt.solve_discrete_problem import solve_MILP_discrete_problem
-from pyomo.contrib.gdpopt.util import _DoNothing
+from pyomo.contrib.gdpopt.util import _DoNothing, is_affine
 from pyomo.core import Block, Constraint, Objective, Var, maximize, value
 from pyomo.gdp import Disjunct
 from pyomo.opt import TerminationCondition as tc
@@ -289,7 +289,7 @@ def init_set_covering(
     # disjunct_list still needs to be covered by the initialization
     disjunct_needs_cover = list(
         any(
-            constr.body.polynomial_degree() not in (0, 1)
+            not is_affine(constr.body)
             for constr in disj.component_data_objects(
                 ctype=Constraint, active=True, descend_into=True
             )

pyomo/contrib/gdpopt/gloa.py

@@ -28,6 +28,7 @@
 from pyomo.contrib.gdpopt.solve_discrete_problem import solve_MILP_discrete_problem
 from pyomo.contrib.gdpopt.util import (
     _add_bigm_constraint_to_transformed_model,
+    is_affine,
     time_code,
 )
 from pyomo.contrib.mcpp.pyomo_mcpp import McCormick as mc, MCPP_Error
@@ -154,7 +155,7 @@ def _add_cuts_to_discrete_problem(
         ):
             disjunctive_var_bounds = disjunctive_bounds(constr.parent_block())
 
-            if constr.body.polynomial_degree() in (1, 0):
+            if is_affine(constr.body):
                 continue
 
             vars_in_constr = list(identify_variables(constr.body))

pyomo/contrib/gdpopt/loa.py

@@ -28,6 +28,7 @@
 from pyomo.contrib.gdpopt.oa_algorithm_utils import _OAAlgorithmMixIn
 from pyomo.contrib.gdpopt.solve_discrete_problem import solve_MILP_discrete_problem
 from pyomo.contrib.gdpopt.util import (
+    is_affine,
     time_code,
     _add_bigm_constraint_to_transformed_model,
 )
@@ -47,7 +48,6 @@
 from pyomo.core.expr.visitor import identify_variables
 from pyomo.gdp import Disjunct
 from pyomo.opt.base import SolverFactory
-from pyomo.repn import generate_standard_repn
 
 MAX_SYMBOLIC_DERIV_SIZE = 1000
 JacInfo = namedtuple('JacInfo', ['mode', 'vars', 'jac'])
@@ -214,7 +214,7 @@ def _add_cuts_to_discrete_problem(
             subproblem_util_block.constraint_list,
         ):
             dual_value = nlp.dual.get(subprob_constr, None)
-            if dual_value is None or generate_standard_repn(constr.body).is_linear():
+            if dual_value is None or is_affine(constr.body):
                 continue
 
             # Determine if the user pre-specified that OA cuts should not be

pyomo/contrib/gdpopt/solve_subproblem.py

@@ -19,8 +19,9 @@
 )
 from pyomo.contrib.gdpopt.util import (
     SuppressInfeasibleWarning,
-    is_feasible,
     get_main_elapsed_time,
+    is_affine,
+    is_feasible,
 )
 from pyomo.core import Constraint, TransformationFactory, Objective, Block
 import pyomo.core.expr as EXPR
@@ -357,7 +358,7 @@ def call_appropriate_subproblem_solver(subprob_util_block, solver, config):
 
     # Is the subproblem linear?
     if not any(
-        constr.body.polynomial_degree() not in (1, 0)
+        not is_affine(constr.body)
         for constr in subprob.component_data_objects(Constraint, active=True)
     ):
         subprob_termination = solve_linear_subproblem(subprob, config, timing)

pyomo/contrib/gdpopt/tests/test_gloa_affine_mutable_param.py

@@ -0,0 +1,247 @@
+# ____________________________________________________________________________________
+#
+# Pyomo: Python Optimization Modeling Objects
+# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC
+# Under the terms of Contract DE-NA0003525 with National Technology and Engineering
+# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this
+# software.  This software is distributed under the 3-clause BSD License.
+# ____________________________________________________________________________________
+
+import logging
+from unittest.mock import patch
+
+import pyomo.environ as pyo
+from pyomo.common import unittest
+from pyomo.contrib.gdpopt.branch_and_bound import GDP_LBB_Solver
+from pyomo.contrib.gdpopt.create_oa_subproblems import add_util_block
+from pyomo.contrib.gdpopt.util import is_affine, move_nonlinear_objective_to_constraints
+from pyomo.gdp import Disjunction
+from pyomo.opt import TerminationCondition
+from pyomo.repn import generate_standard_repn
+
+
+class _StopAfterLBBPreprocessing(Exception):
+    pass
+
+
+def _add_mutable_crf(m):
+    """Add the dimensionless mutable-parameter coefficient from issue #4036."""
+    m.rate = pyo.Param(initialize=0.075, mutable=True)
+    m.years = pyo.Param(initialize=30, mutable=True)
+    m.crf = pyo.Expression(expr=m.rate / (1 - (1 + m.rate) ** (-m.years)))
+
+
+@unittest.skipUnless(
+    pyo.SolverFactory('highs').available(exception_flag=False), 'HiGHS is not available'
+)
+class TestGLOAAffineMutableParameterExpressions(unittest.TestCase):
+    def test_gloa_keeps_affine_objective(self):
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.choose_x = Disjunction(expr=[[m.x == 1], [m.x == 9]])
+        m.obj = pyo.Objective(expr=m.x / m.crf, sense=pyo.maximize)
+
+        self.assertIsNone(m.obj.expr.polynomial_degree())
+        self.assertTrue(generate_standard_repn(m.obj.expr).is_linear())
+        self.assertTrue(is_affine(m.obj.expr))
+
+        result = pyo.SolverFactory('gdpopt.gloa').solve(
+            m, init_algorithm='no_init', iterlim=1, mip_solver='highs'
+        )
+
+        self.assertEqual(
+            result.solver.termination_condition, TerminationCondition.maxIterations
+        )
+        self.assertAlmostEqual(result.problem.upper_bound, 9 / pyo.value(m.crf))
+
+    def test_gloa_does_not_replace_affine_objective(self):
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.obj = pyo.Objective(expr=m.x / m.crf, sense=pyo.maximize)
+        util_block = add_util_block(m)
+        util_block.algebraic_variable_list = []
+
+        original_obj = move_nonlinear_objective_to_constraints(
+            util_block, logging.getLogger(__name__)
+        )
+
+        self.assertIsNone(original_obj)
+        self.assertTrue(m.obj.active)
+        self.assertFalse(hasattr(util_block, 'objective_value'))
+
+    def test_gloa_keeps_affine_constraint(self):
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.limit = pyo.Constraint(expr=m.x / m.crf <= 9 / m.crf)
+        m.choose_x = Disjunction(expr=[[m.x == 1], [m.x == 10]])
+        m.obj = pyo.Objective(expr=m.x, sense=pyo.maximize)
+
+        self.assertIsNone(m.limit.body.polynomial_degree())
+        self.assertTrue(generate_standard_repn(m.limit.body).is_linear())
+        self.assertTrue(is_affine(m.limit.body))
+
+        result = pyo.SolverFactory('gdpopt.gloa').solve(
+            m, init_algorithm='no_init', iterlim=1, mip_solver='highs'
+        )
+
+        self.assertEqual(
+            result.solver.termination_condition, TerminationCondition.maxIterations
+        )
+        self.assertAlmostEqual(result.problem.upper_bound, 1)
+
+    def test_continuous_affine_model_uses_mip_solver(self):
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.limit = pyo.Constraint(expr=m.x / m.crf <= 9 / m.crf)
+        m.obj = pyo.Objective(expr=m.x / m.crf, sense=pyo.maximize)
+
+        result = pyo.SolverFactory('gdpopt.gloa').solve(
+            m, mip_solver='highs', nlp_solver='not_available'
+        )
+
+        self.assertAlmostEqual(result.problem.lower_bound, 9 / pyo.value(m.crf))
+        self.assertAlmostEqual(result.problem.upper_bound, 9 / pyo.value(m.crf))
+        self.assertAlmostEqual(pyo.value(m.x), 9)
+
+    def test_affine_subproblem_uses_mip_solver(self):
+        """Keep an unfixed CRF-coefficient constraint on the LP subproblem."""
+        m = pyo.ConcreteModel()
+        m.selector = pyo.Var(bounds=(0, 1))
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.limit = pyo.Constraint(expr=m.x / m.crf <= 1 / m.crf)
+        m.choose_selector = Disjunction(expr=[[m.selector == 0], [m.selector == 1]])
+        m.obj = pyo.Objective(expr=m.x, sense=pyo.maximize)
+
+        def check_subproblem(solver, subproblem, util_block):
+            self.assertTrue(subproblem.limit.active)
+            self.assertFalse(subproblem.x.fixed)
+            self.assertIsNone(subproblem.limit.body.polynomial_degree())
+            self.assertTrue(is_affine(subproblem.limit.body))
+
+        result = pyo.SolverFactory('gdpopt.gloa').solve(
+            m,
+            init_algorithm='no_init',
+            mip_solver='highs',
+            nlp_solver='not_available',
+            subproblem_presolve=False,
+            call_before_subproblem_solve=check_subproblem,
+        )
+
+        self.assertEqual(
+            result.solver.termination_condition, TerminationCondition.optimal
+        )
+        self.assertAlmostEqual(result.problem.lower_bound, 1)
+        self.assertAlmostEqual(result.problem.upper_bound, 1)
+        self.assertAlmostEqual(pyo.value(m.x), 1)
+
+    def test_set_covering_skips_affine_disjuncts(self):
+        """Do not initialize disjuncts whose constraints are all affine."""
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.choose_x = Disjunction(
+            expr=[[m.x / m.crf <= 1 / m.crf], [m.x / m.crf >= 2 / m.crf]]
+        )
+        m.obj = pyo.Objective(expr=m.x)
+        solver = pyo.SolverFactory('gdpopt.gloa')
+
+        result = solver.solve(
+            m,
+            init_algorithm='set_covering',
+            set_cover_iterlim=1,
+            mip_solver='highs',
+            nlp_solver='not_available',
+        )
+
+        self.assertEqual(solver.initialization_iteration, 0)
+        self.assertEqual(
+            result.solver.termination_condition, TerminationCondition.optimal
+        )
+        self.assertAlmostEqual(pyo.value(m.x), 0)
+
+    def test_gloa_does_not_generate_cuts_for_affine_constraint(self):
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.limit = pyo.Constraint(expr=m.x / m.crf <= 9 / m.crf)
+        m.choose_x = Disjunction(expr=[[m.x == 1], [m.x == 10]])
+        m.obj = pyo.Objective(expr=m.x, sense=pyo.maximize)
+
+        with patch(
+            'pyomo.contrib.gdpopt.gloa.mc',
+            side_effect=AssertionError('affine constraint sent to MC++'),
+        ):
+            result = pyo.SolverFactory('gdpopt.gloa').solve(
+                m, init_algorithm='no_init', mip_solver='highs', nlp_solver='highs'
+            )
+
+        self.assertEqual(
+            result.solver.termination_condition, TerminationCondition.optimal
+        )
+        self.assertAlmostEqual(result.problem.upper_bound, 1)
+
+
+class TestLBBAffineMutableParameterExpressions(unittest.TestCase):
+    def test_lbb_keeps_affine_constraints_in_root_relaxation(self):
+        """Inspect LBB preprocessing without requiring an optional MINLP solver."""
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        _add_mutable_crf(m)
+        m.choose_x = Disjunction(
+            expr=[[m.x / m.crf <= 1 / m.crf], [m.x / m.crf >= 2 / m.crf]]
+        )
+        m.obj = pyo.Objective(expr=m.x)
+
+        def stop_after_preprocessing(solver, node_data, node_model, config):
+            util_block = node_model.component(solver.original_util_block.name)
+            constraints = [
+                constr
+                for disjunct in util_block.disjunct_list
+                for constr in disjunct.component_data_objects(
+                    pyo.Constraint, active=None
+                )
+            ]
+            self.assertEqual(len(constraints), 2)
+            self.assertTrue(all(constr.active for constr in constraints))
+            self.assertEqual(len(util_block.disjunct_to_nonlinear_constraints), 0)
+            raise _StopAfterLBBPreprocessing
+
+        with patch.object(
+            GDP_LBB_Solver, '_prescreen_node', new=stop_after_preprocessing
+        ):
+            with self.assertRaises(_StopAfterLBBPreprocessing):
+                pyo.SolverFactory('gdpopt.lbb').solve(m, minlp_solver='not_available')
+
+
+@unittest.skipUnless(
+    pyo.SolverFactory('glpk').available(exception_flag=False), 'GLPK is not available'
+)
+class TestGLOAAffineExprIf(unittest.TestCase):
+    def test_gloa_keeps_affine_expr_if_constraint(self):
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(bounds=(0, 10))
+        m.sw = pyo.Var(bounds=(0, 1), initialize=1)
+        m.sw.fix(1)
+        m.limit = pyo.Constraint(
+            expr=pyo.Expr_if(IF=m.sw >= 0.5, THEN=m.x, ELSE=2 * m.x) <= 1
+        )
+        m.choose_x = Disjunction(expr=[[m.x >= 0], [m.x == 10]])
+        m.obj = pyo.Objective(expr=m.x, sense=pyo.maximize)
+
+        self.assertEqual(m.limit.body.polynomial_degree(), 1)
+        self.assertFalse(generate_standard_repn(m.limit.body).is_linear())
+        self.assertTrue(is_affine(m.limit.body))
+
+        result = pyo.SolverFactory('gdpopt.gloa').solve(
+            m, init_algorithm='no_init', iterlim=1, mip_solver='glpk'
+        )
+
+        self.assertEqual(
+            result.solver.termination_condition, TerminationCondition.maxIterations
+        )
+        self.assertAlmostEqual(result.problem.upper_bound, 1)