Merge branch 'main' into finalize-release-6.9.2

Author: mrmundt

pyomo/contrib/fbbt/fbbt.py

@@ -1105,7 +1105,9 @@ def _before_external_function(visitor, child):
 def _register_new_before_child_handler(visitor, child):
     handlers = _before_child_handlers
     child_type = child.__class__
-    if child.is_variable_type():
+    if child_type in native_types:
+        handlers[child_type] = _before_constant
+    elif child.is_variable_type():
         handlers[child_type] = _before_var
     elif not child.is_potentially_variable():
         handlers[child_type] = _before_NPV

pyomo/contrib/fbbt/tests/test_fbbt.py

@@ -11,7 +11,7 @@
 
 import pyomo.common.unittest as unittest
 import pyomo.environ as pyo
-from pyomo.contrib.fbbt.fbbt import fbbt, compute_bounds_on_expr
+from pyomo.contrib.fbbt.fbbt import fbbt, compute_bounds_on_expr, _before_child_handlers
 from pyomo.common.dependencies import numpy as np, numpy_available
 from pyomo.common.fileutils import find_library
 from pyomo.common.log import LoggingIntercept
@@ -1363,3 +1363,30 @@ def test_ranged_expression(self):
         self.assertEqual(m.l.bounds, (2, 7))
         self.assertEqual(m.x.bounds, (3, 7))
         self.assertEqual(m.u.bounds, (3, 8))
+
+    @unittest.skipUnless(numpy_available, "Test requires numpy")
+    def test_numpy_leaves(self):
+        m = pyo.ConcreteModel()
+        m.l = pyo.Var(bounds=(2, None))
+        m.x = pyo.Var()
+        m.u = pyo.Var(bounds=(None, 8))
+        m.c = pyo.Constraint(
+            expr=pyo.inequality(m.l + np.int32(1), m.x, m.u - np.float64(1))
+        )
+
+        # Remove the numpy types so we can test that automatic numeric
+        # type registrations
+        old = [(t, _before_child_handlers.pop(t, None)) for t in (np.int32, np.float64)]
+
+        try:
+            self.tightener(m)
+            self.tightener(m)
+            self.assertEqual(m.l.bounds, (2, 6.0))
+            self.assertEqual(m.x.bounds, (3, 7.0))
+            self.assertEqual(m.u.bounds, (4, 8.0))
+        finally:
+            for t, fcn in old:
+                if fcn is None:
+                    _before_child_handlers.pop(t, None)
+                else:
+                    _before_child_handlers[t] = fcn

pyomo/contrib/pyros/master_problem_methods.py

@@ -30,6 +30,7 @@
     call_solver,
     DR_POLISHING_PARAM_PRODUCT_ZERO_TOL,
     enforce_dr_degree,
+    get_all_first_stage_eq_cons,
     get_dr_expression,
     check_time_limit_reached,
     generate_all_decision_rule_var_data_objects,
@@ -137,7 +138,7 @@ def add_scenario_block_to_master_problem(
         new_blk = master_model.scenarios[scenario_idx]
         for con in new_blk.first_stage.inequality_cons.values():
             con.deactivate()
-        for con in new_blk.first_stage.equality_cons.values():
+        for con in get_all_first_stage_eq_cons(new_blk):
             con.deactivate()
 
 
@@ -588,13 +589,21 @@ def minimize_dr_vars(master_data):
 def get_master_dr_degree(master_data):
     """
     Determine DR polynomial degree to enforce based on
-    the iteration number.
+    the iteration number and/or the presence of first-stage
+    equality constraints that depend on the decision rule variables.
 
-    Currently, the degree is set to:
+    If there are first-stage equality constraints that depend
+    on the decision rule variables, such as equalities derived
+    from coefficient matching or discretization of state-variable
+    independent equalities, then the degree is set to
+    ``config.decision_rule_order``.
+
+    Otherwise, the degree is set to:
 
     - 0 if iteration number is 0
     - min(1, config.decision_rule_order) if iteration number
-      otherwise does not exceed number of uncertain parameters
+      otherwise does not exceed number of effective
+      uncertain parameters
     - min(2, config.decision_rule_order) otherwise.
 
     Parameters
@@ -607,9 +616,13 @@ def get_master_dr_degree(master_data):
     int
         DR order, or polynomial degree, to enforce.
     """
+    nom_scenario_blk = master_data.master_model.scenarios[0, 0]
+    if nom_scenario_blk.first_stage.dr_dependent_equality_cons:
+        return master_data.config.decision_rule_order
+
     if master_data.iteration == 0:
         return 0
-    elif master_data.iteration <= len(master_data.config.uncertain_params):
+    elif master_data.iteration <= len(nom_scenario_blk.effective_uncertain_params):
         return min(1, master_data.config.decision_rule_order)
     else:
         return min(2, master_data.config.decision_rule_order)

pyomo/contrib/pyros/separation_problem_methods.py

@@ -38,6 +38,7 @@
     ABS_CON_CHECK_FEAS_TOL,
     call_solver,
     check_time_limit_reached,
+    get_all_first_stage_eq_cons,
 )
 
 
@@ -121,7 +122,7 @@ def construct_separation_problem(model_data):
     # fix/deactivate all nonadjustable components
     for var in separation_model.all_nonadjustable_variables:
         var.fix()
-    for fs_eqcon in separation_model.first_stage.equality_cons.values():
+    for fs_eqcon in get_all_first_stage_eq_cons(separation_model):
         fs_eqcon.deactivate()
     for fs_ineqcon in separation_model.first_stage.inequality_cons.values():
         fs_ineqcon.deactivate()

pyomo/contrib/pyros/tests/test_grcs.py

@@ -1997,6 +1997,58 @@ def test_pyros_certain_params_ipopt_degrees_of_freedom(self):
         )
         self.assertEqual(m.x.value, 1)
 
+    @parameterized.expand([[True, 1], [True, 2], [False, 1], [False, 2]])
+    def test_two_stage_set_nonstatic_dr_robust_opt(self, use_discrete_set, dr_order):
+        """
+        Test problems that are sensitive to the DR order efficiency.
+
+        If the efficiency is not switched off properly, then
+        PyROS may terminate prematurely with a(n inaccurate)
+        robust infeasibility status.
+        """
+        m = ConcreteModel()
+        m.x = Var(bounds=[-2, 2], initialize=0)
+        m.z = Var(bounds=[-10, 10], initialize=0)
+        m.q = Param(initialize=2, mutable=True)
+        m.obj = Objective(expr=m.x + m.z, sense=maximize)
+        # when uncertainty set is discrete, the
+        # preprocessor should write out this constraint for
+        # each scenario as a first-stage constraint
+        # otherwise, coefficient matching constraint
+        # requires only the affine DR coefficient be nonzero
+        m.xz_con = Constraint(expr=m.z == m.q)
+
+        uncertainty_set = (
+            DiscreteScenarioSet([[2], [3]]) if use_discrete_set else BoxSet([[2, 3]])
+        )
+        baron = SolverFactory("baron")
+        res = SolverFactory("pyros").solve(
+            model=m,
+            first_stage_variables=m.x,
+            second_stage_variables=m.z,
+            uncertain_params=m.q,
+            uncertainty_set=uncertainty_set,
+            local_solver=baron,
+            global_solver=baron,
+            solve_master_globally=True,
+            bypass_local_separation=True,
+            decision_rule_order=dr_order,
+            objective_focus="worst_case",
+        )
+        self.assertEqual(
+            # DR efficiency should have been switched off due to
+            # DR-dependent equalities, so robust optimal
+            # if the DR efficiency was not switched off, then
+            # robust infeasibililty would have been prematurely reported
+            res.pyros_termination_condition,
+            pyrosTerminationCondition.robust_optimal,
+        )
+        self.assertEqual(res.iterations, 1)
+        # optimal solution evaluated under worst-case scenario
+        self.assertAlmostEqual(res.final_objective_value, 4, places=4)
+        self.assertAlmostEqual(m.x.value, 2, places=4)
+        self.assertAlmostEqual(m.z.value, 2, places=4)
+
 
 @unittest.skipUnless(baron_available, "BARON not available")
 class TestReformulateSecondStageEqualitiesDiscrete(unittest.TestCase):
@@ -2145,12 +2197,9 @@ def test_two_stage_discrete_set_rank2_affine_dr(self):
             res.pyros_termination_condition, pyrosTerminationCondition.robust_optimal
         )
         self.assertEqual(res.iterations, 1)
-        self.assertAlmostEqual(res.final_objective_value, 0, places=4)
-        # note: this solution is suboptimal (in the context of nonstatic DRs),
-        #       but follows from the current efficiency for DRs
-        #       (i.e. in first iteration, static DRs required)
-        self.assertAlmostEqual(m.x.value, 0, places=4)
-        self.assertAlmostEqual(m.z.value, 0, places=4)
+        self.assertAlmostEqual(res.final_objective_value, 2, places=4)
+        self.assertAlmostEqual(m.x.value, 4, places=4)
+        self.assertAlmostEqual(m.z.value, -2, places=4)
 
     def test_two_stage_discrete_set_fullrank_affine_dr(self):
         m = self.build_two_stage_model()

pyomo/contrib/pyros/tests/test_master.py

@@ -39,6 +39,7 @@
 from pyomo.contrib.pyros.util import (
     ModelData,
     preprocess_model_data,
+    get_all_first_stage_eq_cons,
     ObjectiveType,
     time_code,
     TimingData,
@@ -195,8 +196,8 @@ def test_add_scenario_block_to_master(self):
             )
 
         nadj_eq_con_zip = zip(
-            master_model.scenarios[0, 0].first_stage.equality_cons.values(),
-            master_model.scenarios[0, 1].first_stage.equality_cons.values(),
+            get_all_first_stage_eq_cons(master_model.scenarios[0, 0]),
+            get_all_first_stage_eq_cons(master_model.scenarios[0, 1]),
         )
         for eq_con_00, eq_con_01 in nadj_eq_con_zip:
             self.assertIsNot(