Merge pull request Pyomo#3562 from shermanjasonaf/pyros-discrete-coefficient-matching

Modify PyROS DR Order Efficiency

Author: jsiirola

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(