Revise PyROS second-stage equality reformulation under discrete (scenario-based) uncertainty

pyomo/contrib/pyros/CHANGELOG.txt

@@ -3,6 +3,13 @@ PyROS CHANGELOG
 ===============
 
 
+-------------------------------------------------------------------------------
+PyROS 1.3.7    06 Mar 2025
+-------------------------------------------------------------------------------
+- Modify reformulation of state-variable independent second-stage
+  equality constraints for problems with discrete uncertainty sets
+
+
 -------------------------------------------------------------------------------
 PyROS 1.3.6    06 Mar 2025
 -------------------------------------------------------------------------------

pyomo/contrib/pyros/pyros.py

@@ -33,7 +33,7 @@
 )
 
 
-__version__ = "1.3.6"
+__version__ = "1.3.7"
 
 
 default_pyros_solver_logger = setup_pyros_logger()

pyomo/contrib/pyros/tests/test_grcs.py

@@ -1998,6 +1998,189 @@ def test_pyros_certain_params_ipopt_degrees_of_freedom(self):
         self.assertEqual(m.x.value, 1)
 
 
+@unittest.skipUnless(baron_available, "BARON not available")
+class TestReformulateSecondStageEqualitiesDiscrete(unittest.TestCase):
+    """
+    Test behavior of PyROS solver when the uncertainty set is
+    discrete and there are second-stage
+    equality constraints that are state-variable independent,
+    and therefore, subject to reformulation.
+    """
+
+    def build_single_stage_model(self):
+        m = ConcreteModel()
+        m.x = Var(range(3), bounds=[-2, 2], initialize=0)
+        m.q = Param(range(3), initialize=0, mutable=True)
+        m.c = Param(range(3), initialize={0: 1, 1: 0, 2: 1})
+        m.obj = Objective(expr=sum(m.x[i] * m.c[i] for i in m.x), sense=maximize)
+        # when uncertainty set is discrete, the
+        # preprocessor should write out this constraint for
+        # each scenario as a first-stage constraint
+        m.xq_con = Constraint(expr=sum(m.x[i] * m.q[i] for i in m.x) == 0)
+        return m
+
+    def build_two_stage_model(self):
+        m = ConcreteModel()
+        m.x = Var(bounds=[None, None], initialize=0)
+        m.z = Var(bounds=[-2, 2], 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
+        m.xz_con = Constraint(expr=m.x + m.q * m.z == 0)
+        return m
+
+    def test_single_stage_discrete_set_fullrank(self):
+        m = self.build_single_stage_model()
+        uncertainty_set = DiscreteScenarioSet(
+            # reformulating second-stage equality for these scenarios
+            # should result in first-stage equalities finally being
+            # (full-column-rank matrix) @ (x) == 0
+            # so x=0 is sole robust feasible solution
+            scenarios=[
+                [0] * len(m.q),
+                [1] * len(m.q),
+                list(range(1, len(m.q) + 1)),
+                [(idx + 1) ** 2 for idx in m.q],
+            ]
+        )
+        baron = SolverFactory("baron")
+        res = SolverFactory("pyros").solve(
+            model=m,
+            first_stage_variables=m.x,
+            second_stage_variables=[],
+            uncertain_params=m.q,
+            uncertainty_set=uncertainty_set,
+            local_solver=baron,
+            global_solver=baron,
+            solve_master_globally=True,
+            bypass_local_separation=True,
+            objective_focus="worst_case",
+        )
+        self.assertEqual(
+            res.pyros_termination_condition, pyrosTerminationCondition.robust_optimal
+        )
+        self.assertEqual(res.iterations, 1)
+        self.assertAlmostEqual(res.final_objective_value, 0, places=4)
+        self.assertAlmostEqual(m.x[0].value, 0, places=4)
+        self.assertAlmostEqual(m.x[1].value, 0, places=4)
+        self.assertAlmostEqual(m.x[2].value, 0, places=4)
+
+    def test_single_stage_discrete_set_rank2(self):
+        m = self.build_single_stage_model()
+        uncertainty_set = DiscreteScenarioSet(
+            # reformulating second-stage equality for these scenarios
+            # should make the optimal solution unique
+            scenarios=[[0] * len(m.q), [1] * len(m.q), [(idx + 1) ** 2 for idx in m.q]]
+        )
+        baron = SolverFactory("baron")
+        res = SolverFactory("pyros").solve(
+            model=m,
+            first_stage_variables=m.x,
+            second_stage_variables=[],
+            uncertain_params=m.q,
+            uncertainty_set=uncertainty_set,
+            local_solver=baron,
+            global_solver=baron,
+            solve_master_globally=True,
+            bypass_local_separation=True,
+            objective_focus="worst_case",
+        )
+        self.assertEqual(
+            res.pyros_termination_condition, pyrosTerminationCondition.robust_optimal
+        )
+        self.assertEqual(res.iterations, 1)
+        self.assertAlmostEqual(res.final_objective_value, 2, places=4)
+        # optimal solution is unique
+        self.assertAlmostEqual(m.x[0].value, 5 / 4, places=4)
+        self.assertAlmostEqual(m.x[1].value, -2, places=4)
+        self.assertAlmostEqual(m.x[2].value, 3 / 4, places=4)
+
+    def test_single_stage_discrete_set_rank1(self):
+        m = self.build_single_stage_model()
+        uncertainty_set = DiscreteScenarioSet(
+            scenarios=[[0] * len(m.q), [2] * len(m.q), [3] * len(m.q)]
+        )
+        baron = SolverFactory("baron")
+        res = SolverFactory("pyros").solve(
+            model=m,
+            first_stage_variables=m.x,
+            second_stage_variables=[],
+            uncertain_params=m.q,
+            uncertainty_set=uncertainty_set,
+            local_solver=baron,
+            global_solver=baron,
+            solve_master_globally=True,
+            bypass_local_separation=True,
+            objective_focus="worst_case",
+        )
+        self.assertEqual(
+            res.pyros_termination_condition, pyrosTerminationCondition.robust_optimal
+        )
+        self.assertEqual(res.iterations, 1)
+        self.assertAlmostEqual(res.final_objective_value, 2, places=4)
+        # subject to these scenarios, the optimal solution is non-unique,
+        # but should satisfy this check
+        self.assertAlmostEqual(m.x[1].value, -2, places=4)
+
+    def test_two_stage_discrete_set_rank2_affine_dr(self):
+        m = self.build_two_stage_model()
+        uncertainty_set = DiscreteScenarioSet([[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=1,
+            objective_focus="worst_case",
+        )
+        self.assertEqual(
+            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)
+
+    def test_two_stage_discrete_set_fullrank_affine_dr(self):
+        m = self.build_two_stage_model()
+        uncertainty_set = DiscreteScenarioSet([[2], [3], [4]])
+        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=1,
+            objective_focus="worst_case",
+        )
+        self.assertEqual(
+            res.pyros_termination_condition, pyrosTerminationCondition.robust_optimal
+        )
+        self.assertEqual(res.iterations, 1)
+        self.assertAlmostEqual(res.final_objective_value, 0, places=4)
+        # the second-stage equalities are a full rank linear system
+        # in x and the DR variables, with RHS 0, so all
+        # variables must be 0
+        self.assertAlmostEqual(m.x.value, 0, places=4)
+        self.assertAlmostEqual(m.z.value, 0, places=4)
+
+
 @unittest.skipUnless(ipopt_available, "IPOPT not available.")
 class TestPyROSVarsAsUncertainParams(unittest.TestCase):
     """

pyomo/contrib/pyros/tests/test_preprocessor.py

@@ -36,10 +36,17 @@
     Block,
 )
 from pyomo.core.base.set_types import NonNegativeReals, NonPositiveReals, Reals
-from pyomo.core.expr import LinearExpression, log, sin, exp, RangedExpression
+from pyomo.core.expr import (
+    LinearExpression,
+    log,
+    sin,
+    exp,
+    RangedExpression,
+    SumExpression,
+)
 from pyomo.core.expr.compare import assertExpressionsEqual
 
-from pyomo.contrib.pyros.uncertainty_sets import BoxSet
+from pyomo.contrib.pyros.uncertainty_sets import BoxSet, DiscreteScenarioSet
 from pyomo.contrib.pyros.util import (
     ModelData,
     ObjectiveType,
@@ -1942,13 +1949,13 @@ class TestReformulateStateVarIndependentEqCons(unittest.TestCase):
     state variable-independent second-stage equality constraints.
     """
 
-    def setup_test_model_data(self):
+    def setup_test_model_data(self, uncertainty_set=None):
         """
         Set up simple test model for testing the reformulation
         routine.
         """
         model_data = Bunch()
-        model_data.config = Bunch()
+        model_data.config = Bunch(uncertainty_set=uncertainty_set or BoxSet([[0, 1]]))
         model_data.working_model = working_model = ConcreteModel()
         model_data.working_model.user_model = m = Block()
 
@@ -2155,6 +2162,83 @@ def test_reformulate_nonlinear_state_var_independent_eq_con(self):
             model_data.separation_priority_order["reform_upper_bound_from_eq_con"], 0
         )
 
+    def test_reformulate_equality_cons_discrete_set(self):
+        """
+        Test routine for reformulating state-variable-independent
+        second-stage equality constraints under scenario-based
+        uncertainty works as expected.
+        """
+        model_data = self.setup_test_model_data(
+            uncertainty_set=DiscreteScenarioSet([[0], [0.7]])
+        )
+        model_data.separation_priority_order = dict()
+
+        model_data.config.decision_rule_order = 1
+        model_data.config.progress_logger = logging.getLogger(
+            self.test_reformulate_nonlinear_state_var_independent_eq_con.__name__
+        )
+        model_data.config.progress_logger.setLevel(logging.DEBUG)
+
+        add_decision_rule_variables(model_data)
+        add_decision_rule_constraints(model_data)
+
+        ep = model_data.working_model.effective_var_partitioning
+        model_data.working_model.all_nonadjustable_variables = list(
+            ep.first_stage_variables
+            + list(model_data.working_model.first_stage.decision_rule_var_0.values())
+        )
+
+        wm = model_data.working_model
+        m = model_data.working_model.user_model
+        wm.second_stage.equality_cons["eq_con_2"].set_value(m.u * (m.x1 - 1) == 0)
+
+        robust_infeasible = reformulate_state_var_independent_eq_cons(model_data)
+
+        # check constraint partitioning updated as expected
+        self.assertFalse(robust_infeasible)
+        self.assertFalse(wm.second_stage.equality_cons)
+        self.assertEqual(len(wm.second_stage.inequality_cons), 1)
+        self.assertEqual(len(wm.first_stage.equality_cons), 4)
+
+        self.assertTrue(wm.first_stage.equality_cons["scenario_0_eq_con"].active)
+        self.assertTrue(wm.first_stage.equality_cons["scenario_1_eq_con"].active)
+        self.assertTrue(wm.first_stage.equality_cons["scenario_0_eq_con_2"].active)
+        self.assertTrue(wm.first_stage.equality_cons["scenario_1_eq_con_2"].active)
+
+        # expressions for the new opposing inequalities
+        # and coefficient matching constraint
+        dr_vars = list(wm.first_stage.decision_rule_vars[0].values())
+        assertExpressionsEqual(
+            self,
+            wm.first_stage.equality_cons["scenario_0_eq_con"].expr,
+            (
+                0 * SumExpression([dr_vars[0] + 0 * dr_vars[1], -1])
+                + 0 * (m.x1**3 + 0.5)
+                - ((0 * m.u_cert * m.x1) * (dr_vars[0] + 0 * dr_vars[1]))
+                == (0 * (m.x1 + 2))
+            ),
+        )
+        assertExpressionsEqual(
+            self,
+            wm.first_stage.equality_cons["scenario_1_eq_con"].expr,
+            (
+                (0.7**2) * SumExpression([dr_vars[0] + 0.7 * dr_vars[1], -1])
+                + 0.7 * (m.x1**3 + 0.5)
+                - ((5 * 0.7 * m.u_cert * m.x1) * (dr_vars[0] + 0.7 * dr_vars[1]))
+                == (-0.7 * (m.x1 + 2))
+            ),
+        )
+        assertExpressionsEqual(
+            self,
+            wm.first_stage.equality_cons["scenario_0_eq_con_2"].expr,
+            0 * (m.x1 - 1) == 0,
+        )
+        assertExpressionsEqual(
+            self,
+            wm.first_stage.equality_cons["scenario_1_eq_con_2"].expr,
+            0.7 * (m.x1 - 1) == 0,
+        )
+
     def test_coefficient_matching_robust_infeasible_proof(self):
         """
         Test coefficient matching detects robust infeasibility