diff --git a/pyomo/contrib/pyros/master_problem_methods.py b/pyomo/contrib/pyros/master_problem_methods.py index 69e79577b09..50ebe78ec2b 100644 --- a/pyomo/contrib/pyros/master_problem_methods.py +++ b/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) diff --git a/pyomo/contrib/pyros/separation_problem_methods.py b/pyomo/contrib/pyros/separation_problem_methods.py index ba9881596a8..5d1953a5a51 100644 --- a/pyomo/contrib/pyros/separation_problem_methods.py +++ b/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() diff --git a/pyomo/contrib/pyros/tests/test_grcs.py b/pyomo/contrib/pyros/tests/test_grcs.py index eaf29378161..e4ce53df14c 100644 --- a/pyomo/contrib/pyros/tests/test_grcs.py +++ b/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() diff --git a/pyomo/contrib/pyros/tests/test_master.py b/pyomo/contrib/pyros/tests/test_master.py index 9480cc424d2..4ef4b17157d 100644 --- a/pyomo/contrib/pyros/tests/test_master.py +++ b/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( diff --git a/pyomo/contrib/pyros/tests/test_preprocessor.py b/pyomo/contrib/pyros/tests/test_preprocessor.py index 558cddf4010..8997e2a2707 100644 --- a/pyomo/contrib/pyros/tests/test_preprocessor.py +++ b/pyomo/contrib/pyros/tests/test_preprocessor.py @@ -50,6 +50,7 @@ from pyomo.contrib.pyros.util import ( ModelData, ObjectiveType, + get_all_first_stage_eq_cons, get_effective_var_partitioning, get_var_certain_uncertain_bounds, get_var_bound_pairs, @@ -476,7 +477,8 @@ def test_setup_working_model(self): # constraint partitioning initialization self.assertFalse(working_model.first_stage.inequality_cons) - self.assertFalse(working_model.first_stage.equality_cons) + self.assertFalse(working_model.first_stage.dr_dependent_equality_cons) + self.assertFalse(working_model.first_stage.dr_independent_equality_cons) self.assertFalse(working_model.second_stage.inequality_cons) self.assertFalse(working_model.second_stage.equality_cons) @@ -1291,7 +1293,12 @@ def build_simple_test_model_data(self): model_data.working_model.effective_uncertain_params = [m.q] model_data.working_model.first_stage = Block() - model_data.working_model.first_stage.equality_cons = Constraint(Any) + model_data.working_model.first_stage.dr_dependent_equality_cons = Constraint( + Any + ) + model_data.working_model.first_stage.dr_independent_equality_cons = Constraint( + Any + ) model_data.working_model.second_stage = Block() model_data.working_model.second_stage.equality_cons = Constraint(Any) @@ -1321,7 +1328,7 @@ def test_standardize_equality_constraints(self): standardize_equality_constraints(model_data) - first_stage_eq_cons = working_model.first_stage.equality_cons + first_stage_eq_cons = working_model.first_stage.dr_independent_equality_cons second_stage_eq_cons = working_model.second_stage.equality_cons self.assertEqual(len(first_stage_eq_cons), 2) @@ -1983,7 +1990,8 @@ def setup_test_model_data(self, uncertainty_set=None): working_model.effective_uncertain_params = [m.u] working_model.first_stage = Block() - working_model.first_stage.equality_cons = Constraint(Any) + working_model.first_stage.dr_dependent_equality_cons = Constraint(Any) + working_model.first_stage.dr_independent_equality_cons = Constraint(Any) working_model.second_stage = Block() working_model.second_stage.equality_cons = Constraint(Any) working_model.second_stage.inequality_cons = Constraint(Any) @@ -2038,9 +2046,12 @@ def test_coefficient_matching_correct_constraints_added(self): ), ) - first_stage_eq_cons = model_data.working_model.first_stage.equality_cons + fs_blk = model_data.working_model.first_stage + dr_dependent_fs_eq_cons = fs_blk.dr_dependent_equality_cons + dr_independent_fs_eq_cons = fs_blk.dr_independent_equality_cons + self.assertFalse(dr_dependent_fs_eq_cons) self.assertEqual( - len(first_stage_eq_cons), + len(dr_independent_fs_eq_cons), 3, msg="Number of coefficient matching constraints not as expected.", ) @@ -2050,17 +2061,17 @@ def test_coefficient_matching_correct_constraints_added(self): assertExpressionsEqual( self, - first_stage_eq_cons["coeff_matching_eq_con_coeff_1"].expr, + dr_independent_fs_eq_cons["coeff_matching_eq_con_coeff_1"].expr, m.x1**3 + 0.5 + 5 * m.x1 * m.x2 * (-1) + (-1) * (m.x1 + 2) * (-1) == 0, ) assertExpressionsEqual( self, - first_stage_eq_cons["coeff_matching_eq_con_coeff_2"].expr, + dr_independent_fs_eq_cons["coeff_matching_eq_con_coeff_2"].expr, m.x2 - 1 == 0, ) assertExpressionsEqual( self, - first_stage_eq_cons["coeff_matching_eq_con_2_coeff_1"].expr, + dr_independent_fs_eq_cons["coeff_matching_eq_con_2_coeff_1"].expr, m.x2 - 1 == 0, ) @@ -2112,16 +2123,19 @@ def test_reformulate_nonlinear_state_var_independent_eq_con(self): ), ) + fs_blk = model_data.working_model.first_stage + dr_independent_fs_eq_cons = fs_blk.dr_independent_equality_cons # check constraint partitioning updated as expected self.assertFalse(wm.second_stage.equality_cons) self.assertEqual(len(wm.second_stage.inequality_cons), 3) - self.assertEqual(len(wm.first_stage.equality_cons), 1) + self.assertEqual(len(wm.first_stage.dr_independent_equality_cons), 1) + self.assertFalse(wm.first_stage.dr_dependent_equality_cons) second_stage_ineq_cons = wm.second_stage.inequality_cons self.assertTrue(second_stage_ineq_cons["reform_lower_bound_from_eq_con"].active) self.assertTrue(second_stage_ineq_cons["reform_upper_bound_from_eq_con"].active) self.assertTrue( - wm.first_stage.equality_cons["coeff_matching_eq_con_2_coeff_1"].active + dr_independent_fs_eq_cons["coeff_matching_eq_con_2_coeff_1"].active ) # expressions for the new opposing inequalities @@ -2150,7 +2164,7 @@ def test_reformulate_nonlinear_state_var_independent_eq_con(self): ) assertExpressionsEqual( self, - wm.first_stage.equality_cons["coeff_matching_eq_con_2_coeff_1"].expr, + dr_independent_fs_eq_cons["coeff_matching_eq_con_2_coeff_1"].expr, m.x1 - 1 == 0, ) @@ -2195,22 +2209,25 @@ def test_reformulate_equality_cons_discrete_set(self): robust_infeasible = reformulate_state_var_independent_eq_cons(model_data) # check constraint partitioning updated as expected + dr_dependent_equality_cons = wm.first_stage.dr_dependent_equality_cons + dr_independent_equality_cons = wm.first_stage.dr_independent_equality_cons 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.assertEqual(len(wm.first_stage.dr_dependent_equality_cons), 2) + self.assertEqual(len(wm.first_stage.dr_independent_equality_cons), 2) - 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) + self.assertTrue(dr_dependent_equality_cons["scenario_0_eq_con"].active) + self.assertTrue(dr_dependent_equality_cons["scenario_1_eq_con"].active) + self.assertTrue(dr_independent_equality_cons["scenario_0_eq_con_2"].active) + self.assertTrue(dr_independent_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, + wm.first_stage.dr_dependent_equality_cons["scenario_0_eq_con"].expr, ( 0 * SumExpression([dr_vars[0] + 0 * dr_vars[1], -1]) + 0 * (m.x1**3 + 0.5) @@ -2220,7 +2237,7 @@ def test_reformulate_equality_cons_discrete_set(self): ) assertExpressionsEqual( self, - wm.first_stage.equality_cons["scenario_1_eq_con"].expr, + wm.first_stage.dr_dependent_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) @@ -2230,12 +2247,12 @@ def test_reformulate_equality_cons_discrete_set(self): ) assertExpressionsEqual( self, - wm.first_stage.equality_cons["scenario_0_eq_con_2"].expr, + wm.first_stage.dr_independent_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, + wm.first_stage.dr_independent_equality_cons["scenario_1_eq_con_2"].expr, 0.7 * (m.x1 - 1) == 0, ) @@ -2247,6 +2264,7 @@ def test_coefficient_matching_robust_infeasible_proof(self): # Write the deterministic Pyomo model model_data = self.setup_test_model_data() m = model_data.working_model.user_model + model_data.working_model.first_stage.decision_rule_vars = [] model_data.working_model.second_stage.equality_cons["eq_con"].set_value( expr=m.u * (m.x1**3 + 0.5) - 5 * m.u * m.x1 * m.x2 @@ -2588,8 +2606,12 @@ def test_preprocessor_constraint_partitioning_nonstatic_dr( ), ) self.assertEqual( - list(working_model.first_stage.equality_cons), - ["eq_con_eq2", "eq_con_eq3"] + coeff_matching_con_names, + list(working_model.first_stage.dr_independent_equality_cons), + ["eq_con_eq2", "eq_con_eq3"], + ) + self.assertEqual( + list(working_model.first_stage.dr_dependent_equality_cons), + coeff_matching_con_names, ) self.assertEqual( list(working_model.second_stage.inequality_cons), @@ -2631,7 +2653,7 @@ def test_preprocessor_constraint_partitioning_nonstatic_dr( ) # verify the constraints are active - for fs_eq_con in working_model.first_stage.equality_cons.values(): + for fs_eq_con in get_all_first_stage_eq_cons(working_model): self.assertTrue(fs_eq_con.active, msg=f"{fs_eq_con.name} inactive") for fs_ineq_con in working_model.first_stage.inequality_cons.values(): self.assertTrue(fs_ineq_con.active, msg=f"{fs_ineq_con.name} inactive") @@ -2713,11 +2735,11 @@ def test_preprocessor_constraint_partitioning_nonstatic_dr( ) assertExpressionsEqual( - self, fs.equality_cons["eq_con_eq2"].expr, m.x1 - m.z1 == 0 + self, fs.dr_independent_equality_cons["eq_con_eq2"].expr, m.x1 - m.z1 == 0 ) assertExpressionsEqual( self, - fs.equality_cons["eq_con_eq3"].expr, + fs.dr_independent_equality_cons["eq_con_eq3"].expr, m.x1**2 + m.x2 * m.q_cert + m.p * m.z2 == m.p, ) if dr_order < 2: @@ -2770,40 +2792,46 @@ def test_preprocessor_coefficient_matching( working_model = model_data.working_model m = model_data.working_model.user_model fs = working_model.first_stage - fs_eqs = working_model.first_stage.equality_cons + dr_dependent_fs_eqs = working_model.first_stage.dr_dependent_equality_cons ss_ineqs = working_model.second_stage.inequality_cons if config.decision_rule_order == 1: # check the constraint expressions of eq1 and z5 bound assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_0"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_0" + ].expr, fs.decision_rule_vars[1][0] == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_1"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_1" + ].expr, fs.decision_rule_vars[1][1] - 1 == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_2"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_2" + ].expr, fs.decision_rule_vars[1][2] == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_eq_con_eq1_coeff_1"].expr, + dr_dependent_fs_eqs["coeff_matching_eq_con_eq1_coeff_1"].expr, # note: the certain parameter was eliminated by # the expression parser fs.decision_rule_vars[0][0] + m.x2 == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_eq_con_eq1_coeff_2"].expr, + dr_dependent_fs_eqs["coeff_matching_eq_con_eq1_coeff_2"].expr, fs.decision_rule_vars[0][1] == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_eq_con_eq1_coeff_3"].expr, + dr_dependent_fs_eqs["coeff_matching_eq_con_eq1_coeff_3"].expr, fs.decision_rule_vars[0][2] == 0, ) if config.decision_rule_order == 2: @@ -2822,32 +2850,44 @@ def test_preprocessor_coefficient_matching( # check coefficient matching constraint expressions assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_0"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_0" + ].expr, fs.decision_rule_vars[1][0] == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_1"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_1" + ].expr, fs.decision_rule_vars[1][1] - 1 == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_2"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_2" + ].expr, fs.decision_rule_vars[1][2] == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_3"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_3" + ].expr, fs.decision_rule_vars[1][3] == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_4"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_4" + ].expr, fs.decision_rule_vars[1][4] == 0, ) assertExpressionsEqual( self, - fs_eqs["coeff_matching_var_z5_uncertain_eq_bound_con_coeff_5"].expr, + dr_dependent_fs_eqs[ + "coeff_matching_var_z5_uncertain_eq_bound_con_coeff_5" + ].expr, fs.decision_rule_vars[1][5] == 0, ) diff --git a/pyomo/contrib/pyros/util.py b/pyomo/contrib/pyros/util.py index 5d9de4f1bf0..257d2529bf5 100644 --- a/pyomo/contrib/pyros/util.py +++ b/pyomo/contrib/pyros/util.py @@ -1823,7 +1823,8 @@ def setup_working_model(model_data, user_var_partitioning): # stagewise blocks for containing stagewise constraints working_model.first_stage = Block() - working_model.first_stage.equality_cons = Constraint(Any) + working_model.first_stage.dr_independent_equality_cons = Constraint(Any) + working_model.first_stage.dr_dependent_equality_cons = Constraint(Any) working_model.first_stage.inequality_cons = Constraint(Any) working_model.second_stage = Block() working_model.second_stage.equality_cons = Constraint(Any) @@ -2028,7 +2029,9 @@ def standardize_equality_constraints(model_data): con.expr ) else: - working_model.first_stage.equality_cons[f"eq_con_{con_rel_name}"] = con.expr + working_model.first_stage.dr_independent_equality_cons[ + f"eq_con_{con_rel_name}" + ] = con.expr # definitely don't want active duplicate con.deactivate() @@ -2238,6 +2241,12 @@ def get_all_nonadjustable_variables(working_model): return [epigraph_var] + decision_rule_vars + effective_first_stage_vars +def get_all_first_stage_eq_cons(working_model): + return list(working_model.first_stage.dr_dependent_equality_cons.values()) + list( + working_model.first_stage.dr_independent_equality_cons.values() + ) + + def get_all_adjustable_variables(working_model): """ Get all variables considered adjustable. @@ -2344,7 +2353,12 @@ def check_time_limit_reached(timing_data, config): def _reformulate_eq_con_scenario_uncertainty( - model_data, discrete_set, ss_eq_con, ss_eq_con_index, ss_var_id_to_dr_expr_map + model_data, + discrete_set, + ss_eq_con, + ss_eq_con_index, + ss_var_id_to_dr_expr_map, + all_dr_vars_set, ): """ Reformulate a second-stage equality constraint that @@ -2373,21 +2387,36 @@ def _reformulate_eq_con_scenario_uncertainty( ss_var_id_to_dr_expr_map : dict Mapping from object IDs of second-stage variables to corresponding decision rule expressions. + all_dr_vars_set : ComponentSet + A set of all the decision rule variables declared on + ``model_data.working_model``. """ working_model = model_data.working_model con_expr_after_dr_substitution = replace_expressions( expr=ss_eq_con.expr, substitution_map=ss_var_id_to_dr_expr_map ) + vars_in_coeff_expr = ComponentSet( + identify_variables(con_expr_after_dr_substitution) + ) + has_con_dr_vars = vars_in_coeff_expr & all_dr_vars_set + indexed_con_to_update = ( + working_model.first_stage.dr_dependent_equality_cons + if has_con_dr_vars + else working_model.first_stage.dr_independent_equality_cons + ) + scenarios_enum = enumerate(discrete_set.scenarios) for sc_idx, scenario in scenarios_enum: - working_model.first_stage.equality_cons[ - f"scenario_{sc_idx}_{ss_eq_con_index}" - ] = replace_expressions( - expr=con_expr_after_dr_substitution, - substitution_map={ - id(param): scenario_val - for param, scenario_val in zip(working_model.uncertain_params, scenario) - }, + indexed_con_to_update[f"scenario_{sc_idx}_{ss_eq_con_index}"] = ( + replace_expressions( + expr=con_expr_after_dr_substitution, + substitution_map={ + id(param): scenario_val + for param, scenario_val in zip( + working_model.uncertain_params, scenario + ) + }, + ) ) del working_model.second_stage.equality_cons[ss_eq_con_index] @@ -2400,6 +2429,7 @@ def _reformulate_eq_con_continuous_uncertainty( ss_var_id_to_dr_expr_map, uncertain_param_id_to_temp_var_map, originally_unfixed_vars, + all_dr_vars_set, ): """ Reformulate a second-stage equality constraint that @@ -2440,6 +2470,9 @@ def _reformulate_eq_con_continuous_uncertainty( originally_unfixed_vars : list/ComponentSet of VarData Variables of the working model that were originally unfixed. + all_dr_vars_set : ComponentSet + A set of all the decision rule variables declared on + ``model_data.working_model``. Returns ------- @@ -2544,15 +2577,26 @@ def _reformulate_eq_con_continuous_uncertainty( break else: - # coefficient is dependent on model first-stage - # and DR variables. add matching constraint + # coefficient is variable-dependent. + # add matching constraint new_con_name = f"coeff_matching_{ss_eq_con_index}_coeff_{coeff_idx}" - working_model.first_stage.equality_cons[new_con_name] = coeff_expr == 0 - new_con = working_model.first_stage.equality_cons[new_con_name] + vars_in_coeff_expr = ComponentSet(identify_variables(coeff_expr)) + has_expr_dr_vars = vars_in_coeff_expr & all_dr_vars_set + indexed_con_to_update = ( + working_model.first_stage.dr_dependent_equality_cons + if has_expr_dr_vars + else working_model.first_stage.dr_independent_equality_cons + ) + indexed_con_to_update[new_con_name] = coeff_expr == 0 + new_con = indexed_con_to_update[new_con_name] working_model.first_stage.coefficient_matching_cons.append(new_con) + dr_dependence_qual = ( + f"DR variable-{'in' if has_expr_dr_vars else ''}dependent" + ) config.progress_logger.debug( - f"Derived from constraint {ss_eq_con.name!r} a coefficient " + f"Derived from constraint {ss_eq_con.name!r} a " + f"{dr_dependence_qual} coefficient " f"matching constraint named {new_con_name!r} " "with expression: \n " f"{new_con.expr}." @@ -2605,6 +2649,9 @@ def reformulate_state_var_independent_eq_cons(model_data): effective_second_stage_var_set = ComponentSet(ep.second_stage_variables) effective_state_var_set = ComponentSet(ep.state_variables) all_vars_set = ComponentSet(working_model.all_variables) + all_dr_vars_set = ComponentSet( + generate_all_decision_rule_var_data_objects(working_model) + ) originally_unfixed_vars = [var for var in all_vars_set if not var.fixed] # we will need this to substitute DR expressions for @@ -2660,6 +2707,7 @@ def reformulate_state_var_independent_eq_cons(model_data): ss_eq_con_index=con_idx, discrete_set=config.uncertainty_set, ss_var_id_to_dr_expr_map=ss_var_id_to_dr_expr_map, + all_dr_vars_set=all_dr_vars_set, ) else: robust_infeasible = _reformulate_eq_con_continuous_uncertainty( @@ -2672,6 +2720,7 @@ def reformulate_state_var_independent_eq_cons(model_data): uncertain_param_id_to_temp_var_map ), originally_unfixed_vars=originally_unfixed_vars, + all_dr_vars_set=all_dr_vars_set, ) if robust_infeasible: break @@ -2818,11 +2867,12 @@ def log_model_statistics(model_data): # # equality constraints num_eq_cons = ( - len(working_model.first_stage.equality_cons) + len(working_model.first_stage.dr_dependent_equality_cons) + + len(working_model.first_stage.dr_independent_equality_cons) + len(working_model.second_stage.equality_cons) + len(working_model.second_stage.decision_rule_eqns) ) - num_first_stage_eq_cons = len(working_model.first_stage.equality_cons) + num_first_stage_eq_cons = len(get_all_first_stage_eq_cons(working_model)) num_coeff_matching_cons = len(working_model.first_stage.coefficient_matching_cons) num_other_first_stage_eqns = num_first_stage_eq_cons - num_coeff_matching_cons num_second_stage_eq_cons = len(working_model.second_stage.equality_cons)