diff --git a/doc/OnlineDocs/contributed_packages/pyros.rst b/doc/OnlineDocs/contributed_packages/pyros.rst index 133258fb9b8..3ff1bfccf0e 100644 --- a/doc/OnlineDocs/contributed_packages/pyros.rst +++ b/doc/OnlineDocs/contributed_packages/pyros.rst @@ -723,11 +723,11 @@ For this example, we obtain the following price of robustness results: +==========================================+==============================+=============================+ | 0.00 | 35,837,659.18 | 0.00 % | +------------------------------------------+------------------------------+-----------------------------+ - | 0.10 | 36,135,191.59 | 0.82 % | + | 0.10 | 36,135,182.66 | 0.83 % | +------------------------------------------+------------------------------+-----------------------------+ - | 0.20 | 36,437,979.81 | 1.64 % | + | 0.20 | 36,437,979.81 | 1.68 % | +------------------------------------------+------------------------------+-----------------------------+ - | 0.30 | 43,478,190.92 | 17.57 % | + | 0.30 | 43,478,190.91 | 21.32 % | +------------------------------------------+------------------------------+-----------------------------+ | 0.40 | ``robust_infeasible`` | :math:`\text{-----}` | +------------------------------------------+------------------------------+-----------------------------+ @@ -854,10 +854,10 @@ Observe that the log contains the following information: :linenos: ============================================================================== - PyROS: The Pyomo Robust Optimization Solver, v1.2.8. + PyROS: The Pyomo Robust Optimization Solver, v1.2.9. Pyomo version: 6.7.0 Commit hash: unknown - Invoked at UTC 2023-11-03T04:27:42.954101 + Invoked at UTC 2023-12-16T00:00:00.000000 Developed by: Natalie M. Isenberg (1), Jason A. F. Sherman (1), John D. Siirola (2), Chrysanthos E. Gounaris (1) @@ -914,14 +914,11 @@ Observe that the log contains the following information: ------------------------------------------------------------------------------ Itn Objective 1-Stg Shift 2-Stg Shift #CViol Max Viol Wall Time (s) ------------------------------------------------------------------------------ - 0 3.5838e+07 - - 5 1.8832e+04 1.555 - 1 3.5838e+07 2.2045e-12 2.7854e-12 7 3.7766e+04 2.991 - 2 3.6116e+07 1.2324e-01 3.9256e-01 8 1.3466e+06 4.881 - 3 3.6285e+07 5.1968e-01 4.5604e-01 6 4.8734e+03 6.908 - 4 3.6285e+07 2.6524e-13 1.3909e-13 1 3.5036e+01 9.176 - 5 3.6285e+07 2.0167e-13 5.4357e-02 6 2.9103e+00 11.457 - 6 3.6285e+07 2.2335e-12 1.2150e-01 5 4.1726e-01 13.868 - 7 3.6285e+07 2.2340e-12 1.1422e-01 0 9.3279e-10g 20.917 + 0 3.5838e+07 - - 5 1.8832e+04 1.741 + 1 3.5838e+07 3.5184e-15 3.9404e-15 10 4.2516e+06 3.766 + 2 3.5993e+07 1.8105e-01 7.1406e-01 13 5.2004e+06 6.288 + 3 3.6285e+07 5.1968e-01 7.7753e-01 4 1.7892e+04 8.247 + 4 3.6285e+07 9.1166e-13 1.9702e-15 0 7.1157e-10g 11.456 ------------------------------------------------------------------------------ Robust optimal solution identified. ------------------------------------------------------------------------------ @@ -929,22 +926,22 @@ Observe that the log contains the following information: Identifier ncalls cumtime percall % ----------------------------------------------------------- - main 1 20.918 20.918 100.0 + main 1 11.457 11.457 100.0 ------------------------------------------------------ - dr_polishing 7 1.472 0.210 7.0 - global_separation 47 1.239 0.026 5.9 - local_separation 376 9.244 0.025 44.2 - master 8 5.259 0.657 25.1 - master_feasibility 7 0.486 0.069 2.3 - preprocessing 1 0.403 0.403 1.9 - other n/a 2.815 n/a 13.5 + dr_polishing 4 0.682 0.171 6.0 + global_separation 47 1.109 0.024 9.7 + local_separation 235 5.810 0.025 50.7 + master 5 1.353 0.271 11.8 + master_feasibility 4 0.247 0.062 2.2 + preprocessing 1 0.429 0.429 3.7 + other n/a 1.828 n/a 16.0 ====================================================== =========================================================== ------------------------------------------------------------------------------ Termination stats: - Iterations : 8 - Solve time (wall s) : 20.918 + Iterations : 5 + Solve time (wall s) : 11.457 Final objective value : 3.6285e+07 Termination condition : pyrosTerminationCondition.robust_optimal ------------------------------------------------------------------------------ diff --git a/pyomo/contrib/pyros/CHANGELOG.txt b/pyomo/contrib/pyros/CHANGELOG.txt index 7977c37fb95..7d4678f0ba3 100644 --- a/pyomo/contrib/pyros/CHANGELOG.txt +++ b/pyomo/contrib/pyros/CHANGELOG.txt @@ -2,6 +2,18 @@ PyROS CHANGELOG =============== +------------------------------------------------------------------------------- +PyROS 1.2.9 15 Dec 2023 +------------------------------------------------------------------------------- +- Fix DR polishing optimality constraint for case of nominal objective focus +- Use previous separation solution to initialize second-stage and state + variables of new master block; simplify the master feasibility problem +- Use best known solution from master to initialize separation problems + per performance constraint +- Refactor DR variable and constraint declaration routines. +- Refactor DR polishing routine; initialize auxiliary variables + to values they are meant to represent + ------------------------------------------------------------------------------- PyROS 1.2.8 12 Oct 2023 ------------------------------------------------------------------------------- diff --git a/pyomo/contrib/pyros/master_problem_methods.py b/pyomo/contrib/pyros/master_problem_methods.py index 58583c65fba..5477fcc5048 100644 --- a/pyomo/contrib/pyros/master_problem_methods.py +++ b/pyomo/contrib/pyros/master_problem_methods.py @@ -36,7 +36,7 @@ from pyomo.common.modeling import unique_component_name from pyomo.common.timing import TicTocTimer -from pyomo.contrib.pyros.util import TIC_TOC_SOLVE_TIME_ATTR +from pyomo.contrib.pyros.util import TIC_TOC_SOLVE_TIME_ATTR, enforce_dr_degree def initial_construct_master(model_data): @@ -102,6 +102,11 @@ def construct_master_feasibility_problem(model_data, config): Slack variable model. """ + # clone master model. current state: + # - variables for all but newest block are set to values from + # master solution from previous iteration + # - variables for newest block are set to values from separation + # solution chosen in previous iteration model = model_data.master_model.clone() # obtain mapping from master problem to master feasibility @@ -123,53 +128,28 @@ def construct_master_feasibility_problem(model_data, config): obj.deactivate() iteration = model_data.iteration - # first stage vars are already initialized appropriately. - # initialize second-stage DOF variables using DR equation expressions - if model.scenarios[iteration, 0].util.second_stage_variables: - for blk in model.scenarios[iteration, :]: - for eq in blk.util.decision_rule_eqns: - vars_in_dr_eq = ComponentSet(identify_variables(eq.body)) - ssv_set = ComponentSet(blk.util.second_stage_variables) - - # get second-stage var in DR eqn. should only be one var - ssv_in_dr_eq = [var for var in vars_in_dr_eq if var in ssv_set][0] - - # update var value for initialization - # fine since DR eqns are f(d) - z == 0 (not z - f(d) == 0) - ssv_in_dr_eq.set_value(0) - ssv_in_dr_eq.set_value(value(eq.body)) - - # initialize state vars to previous master solution values - if iteration != 0: - stvar_map = get_state_vars(model, [iteration, iteration - 1]) - for current, prev in zip(stvar_map[iteration], stvar_map[iteration - 1]): - current.set_value(value(prev)) - - # constraints to which slacks should be added - # (all the constraints for the current iteration, except the DR eqns) + # add slacks only to inequality constraints for the newest + # master block. these should be the only constraints which + # may have been violated by the previous master and separation + # solution(s) targets = [] for blk in model.scenarios[iteration, :]: - if blk.util.second_stage_variables: - dr_eqs = blk.util.decision_rule_eqns - else: - dr_eqs = list() - targets.extend( [ con for con in blk.component_data_objects( Constraint, active=True, descend_into=True ) - if con not in dr_eqs + if not con.equality ] ) - # retain original constraint exprs (for slack initialization and scaling) + # retain original constraint expressions + # (for slack initialization and scaling) pre_slack_con_exprs = ComponentMap((con, con.body - con.upper) for con in targets) # add slack variables and objective - # inequalities g(v) <= b become g(v) -- s^-<= b - # equalities h(v) == b become h(v) -- s^- + s^+ == b + # inequalities g(v) <= b become g(v) - s^- <= b TransformationFactory("core.add_slack_variables").apply_to(model, targets=targets) slack_vars = ComponentSet( model._core_add_slack_variables.component_data_objects(Var, descend_into=True) @@ -177,8 +157,8 @@ def construct_master_feasibility_problem(model_data, config): # initialize and scale slack variables for con in pre_slack_con_exprs: - # obtain slack vars in updated constraints - # and their coefficients (+/-1) in the constraint expression + # get mapping from slack variables to their (linear) + # coefficients (+/-1) in the updated constraint expressions repn = generate_standard_repn(con.body) slack_var_coef_map = ComponentMap() for idx in range(len(repn.linear_vars)): @@ -187,19 +167,19 @@ def construct_master_feasibility_problem(model_data, config): slack_var_coef_map[var] = repn.linear_coefs[idx] slack_substitution_map = dict() - for slack_var in slack_var_coef_map: - # coefficient determines whether the slack is a +ve or -ve slack + # coefficient determines whether the slack + # is a +ve or -ve slack if slack_var_coef_map[slack_var] == -1: con_slack = max(0, value(pre_slack_con_exprs[con])) else: con_slack = max(0, -value(pre_slack_con_exprs[con])) - # initialize slack var, evaluate scaling coefficient - scaling_coeff = 1 + # initialize slack variable, evaluate scaling coefficient slack_var.set_value(con_slack) + scaling_coeff = 1 - # update expression replacement map + # update expression replacement map for slack scaling slack_substitution_map[id(slack_var)] = scaling_coeff * slack_var # finally, scale slack(s) @@ -307,11 +287,10 @@ def solve_master_feasibility_problem(model_data, config): return results -def minimize_dr_vars(model_data, config): +def construct_dr_polishing_problem(model_data, config): """ - Polish the PyROS decision rule determined for the most - recently solved master problem by minimizing the collective - L1 norm of the vector of all decision rule variables. + Construct DR polishing problem from most recently added + master problem. Parameters ---------- @@ -322,174 +301,153 @@ def minimize_dr_vars(model_data, config): Returns ------- - results : SolverResults - Subordinate solver results for the polishing problem. - polishing_successful : bool - True if polishing model was solved to acceptable level, - False otherwise. + polishing_model : ConcreteModel + Polishing model. + + Note + ---- + Polishing problem is to minimize the L1-norm of the vector of + all decision rule polynomial terms, subject to the original + master problem constraints, with all first-stage variables + (including epigraph) fixed. Optimality of the polished + DR with respect to the master objective is also enforced. """ - # config.progress_logger.info("Executing decision rule variable polishing solve.") - model = model_data.master_model - polishing_model = model.clone() + # clone master problem + master_model = model_data.master_model + polishing_model = master_model.clone() + nominal_polishing_block = polishing_model.scenarios[0, 0] + + # fix first-stage variables (including epigraph, where applicable) + decision_rule_var_set = ComponentSet( + var + for indexed_dr_var in nominal_polishing_block.util.decision_rule_vars + for var in indexed_dr_var.values() + ) + first_stage_vars = nominal_polishing_block.util.first_stage_variables + for var in first_stage_vars: + if var not in decision_rule_var_set: + var.fix() - first_stage_variables = polishing_model.scenarios[0, 0].util.first_stage_variables - decision_rule_vars = polishing_model.scenarios[0, 0].util.decision_rule_vars + # ensure master optimality constraint enforced + if config.objective_focus == ObjectiveType.worst_case: + polishing_model.zeta.fix() + else: + optimal_master_obj_value = value(polishing_model.obj) + polishing_model.nominal_optimality_con = Constraint( + expr=( + nominal_polishing_block.first_stage_objective + + nominal_polishing_block.second_stage_objective + <= optimal_master_obj_value + ) + ) + # deactivate master problem objective polishing_model.obj.deactivate() - index_set = decision_rule_vars[0].index_set() - polishing_model.tau_vars = [] - # ========== - for idx in range(len(decision_rule_vars)): - polishing_model.scenarios[0, 0].add_component( - "polishing_var_" + str(idx), - Var(index_set, initialize=1e6, domain=NonNegativeReals), - ) - polishing_model.tau_vars.append( - getattr(polishing_model.scenarios[0, 0], "polishing_var_" + str(idx)) + + decision_rule_vars = nominal_polishing_block.util.decision_rule_vars + nominal_polishing_block.util.polishing_vars = polishing_vars = [] + for idx, indexed_dr_var in enumerate(decision_rule_vars): + # declare auxiliary 'polishing' variables. + # these are meant to represent the absolute values + # of the terms of DR polynomial + indexed_polishing_var = Var( + list(indexed_dr_var.keys()), domain=NonNegativeReals ) - # ========== - this_iter = polishing_model.scenarios[max(polishing_model.scenarios.keys())[0], 0] - nom_block = polishing_model.scenarios[0, 0] - if config.objective_focus == ObjectiveType.nominal: - obj_val = value( - nom_block.second_stage_objective + nom_block.first_stage_objective + nominal_polishing_block.add_component( + unique_component_name(nominal_polishing_block, f"dr_polishing_var_{idx}"), + indexed_polishing_var, ) - polishing_model.scenarios[0, 0].polishing_constraint = Constraint( - expr=obj_val - >= nom_block.second_stage_objective + nom_block.first_stage_objective + polishing_vars.append(indexed_polishing_var) + + dr_eq_var_zip = zip( + nominal_polishing_block.util.decision_rule_eqns, + polishing_vars, + nominal_polishing_block.util.second_stage_variables, + ) + nominal_polishing_block.util.polishing_abs_val_lb_cons = all_lb_cons = [] + nominal_polishing_block.util.polishing_abs_val_ub_cons = all_ub_cons = [] + for idx, (dr_eq, indexed_polishing_var, ss_var) in enumerate(dr_eq_var_zip): + # set up absolute value constraint components + polishing_absolute_value_lb_cons = Constraint(indexed_polishing_var.index_set()) + polishing_absolute_value_ub_cons = Constraint(indexed_polishing_var.index_set()) + + # add constraints to polishing model + nominal_polishing_block.add_component( + unique_component_name(polishing_model, f"polishing_abs_val_lb_con_{idx}"), + polishing_absolute_value_lb_cons, ) - elif config.objective_focus == ObjectiveType.worst_case: - polishing_model.zeta.fix() # Searching equivalent optimal solutions given optimal zeta - - # === Make absolute value constraints on polishing_vars - polishing_model.scenarios[ - 0, 0 - ].util.absolute_var_constraints = cons = ConstraintList() - uncertain_params = nom_block.util.uncertain_params - if config.decision_rule_order == 1: - for i, tau in enumerate(polishing_model.tau_vars): - for j in range(len(this_iter.util.decision_rule_vars[i])): - if j == 0: - cons.add(-tau[j] <= this_iter.util.decision_rule_vars[i][j]) - cons.add(this_iter.util.decision_rule_vars[i][j] <= tau[j]) - else: - cons.add( - -tau[j] - <= this_iter.util.decision_rule_vars[i][j] - * uncertain_params[j - 1] - ) - cons.add( - this_iter.util.decision_rule_vars[i][j] - * uncertain_params[j - 1] - <= tau[j] - ) - elif config.decision_rule_order == 2: - l = list(range(len(uncertain_params))) - index_pairs = list(it.combinations(l, 2)) - for i, tau in enumerate(polishing_model.tau_vars): - Z = this_iter.util.decision_rule_vars[i] - indices = list(k for k in range(len(Z))) - for r in indices: - if r == 0: - cons.add(-tau[r] <= Z[r]) - cons.add(Z[r] <= tau[r]) - elif r <= len(uncertain_params) and r > 0: - cons.add(-tau[r] <= Z[r] * uncertain_params[r - 1]) - cons.add(Z[r] * uncertain_params[r - 1] <= tau[r]) - elif r <= len(indices) - len(uncertain_params) - 1 and r > len( - uncertain_params - ): - cons.add( - -tau[r] - <= Z[r] - * uncertain_params[ - index_pairs[r - len(uncertain_params) - 1][0] - ] - * uncertain_params[ - index_pairs[r - len(uncertain_params) - 1][1] - ] - ) - cons.add( - Z[r] - * uncertain_params[ - index_pairs[r - len(uncertain_params) - 1][0] - ] - * uncertain_params[ - index_pairs[r - len(uncertain_params) - 1][1] - ] - <= tau[r] - ) - elif r > len(indices) - len(uncertain_params) - 1: - cons.add( - -tau[r] - <= Z[r] - * uncertain_params[ - r - len(index_pairs) - len(uncertain_params) - 1 - ] - ** 2 - ) - cons.add( - Z[r] - * uncertain_params[ - r - len(index_pairs) - len(uncertain_params) - 1 - ] - ** 2 - <= tau[r] - ) - else: - raise NotImplementedError( - "Decision rule variable polishing has not been generalized to decision_rule_order " - + str(config.decision_rule_order) - + "." + nominal_polishing_block.add_component( + unique_component_name(polishing_model, f"polishing_abs_val_ub_con_{idx}"), + polishing_absolute_value_ub_cons, ) - polishing_model.scenarios[0, 0].polishing_obj = Objective( - expr=sum( - sum(tau[j] for j in tau.index_set()) for tau in polishing_model.tau_vars - ) + # update list of absolute value cons + all_lb_cons.append(polishing_absolute_value_lb_cons) + all_ub_cons.append(polishing_absolute_value_ub_cons) + + # get monomials; ensure second-stage variable term excluded + dr_expr_terms = dr_eq.body.args[:-1] + + for dr_eq_term in dr_expr_terms: + dr_var_in_term = dr_eq_term.args[-1] + dr_var_in_term_idx = dr_var_in_term.index() + + # get corresponding polishing variable + polishing_var = indexed_polishing_var[dr_var_in_term_idx] + + # add polishing constraints + polishing_absolute_value_lb_cons[dr_var_in_term_idx] = ( + -polishing_var - dr_eq_term <= 0 + ) + polishing_absolute_value_ub_cons[dr_var_in_term_idx] = ( + dr_eq_term - polishing_var <= 0 + ) + + # if DR var is fixed, then fix corresponding polishing + # variable, and deactivate the absolute value constraints + if dr_var_in_term.fixed: + polishing_var.fix() + polishing_absolute_value_lb_cons[dr_var_in_term_idx].deactivate() + polishing_absolute_value_ub_cons[dr_var_in_term_idx].deactivate() + + # initialize polishing variable to absolute value of + # the DR term. polishing constraints should now be + # satisfied (to equality) at the initial point + polishing_var.set_value(abs(value(dr_eq_term))) + + # polishing problem objective is taken to be 1-norm + # of DR monomials, or equivalently, sum of the polishing + # variables. + polishing_model.polishing_obj = Objective( + expr=sum(sum(polishing_var.values()) for polishing_var in polishing_vars) ) - # === Fix design - for d in first_stage_variables: - d.fix() - - # === Unfix DR vars - num_dr_vars = len( - model.scenarios[0, 0].util.decision_rule_vars[0] - ) # there is at least one dr var - num_uncertain_params = len(config.uncertain_params) - - if model.const_efficiency_applied: - for d in decision_rule_vars: - for i in range(1, num_dr_vars): - d[i].fix(0) - d[0].unfix() - elif model.linear_efficiency_applied: - for d in decision_rule_vars: - d.unfix() - for i in range(num_uncertain_params + 1, num_dr_vars): - d[i].fix(0) - else: - for d in decision_rule_vars: - d.unfix() - - # === Unfix all control var values - for block in polishing_model.scenarios.values(): - for c in block.util.second_stage_variables: - c.unfix() - if model.const_efficiency_applied: - for d in block.util.decision_rule_vars: - for i in range(1, num_dr_vars): - d[i].fix(0) - d[0].unfix() - elif model.linear_efficiency_applied: - for d in block.util.decision_rule_vars: - d.unfix() - for i in range(num_uncertain_params + 1, num_dr_vars): - d[i].fix(0) - else: - for d in block.util.decision_rule_vars: - d.unfix() + return polishing_model + + +def minimize_dr_vars(model_data, config): + """ + Polish decision rule of most recent master problem solution. + + Parameters + ---------- + model_data : MasterProblemData + Master problem data. + config : ConfigDict + PyROS solver settings. + + Returns + ------- + results : SolverResults + Subordinate solver results for the polishing problem. + polishing_successful : bool + True if polishing model was solved to acceptable level, + False otherwise. + """ + # create polishing NLP + polishing_model = construct_dr_polishing_problem( + model_data=model_data, config=config + ) if config.solve_master_globally: solver = config.global_solver @@ -498,11 +456,10 @@ def minimize_dr_vars(model_data, config): config.progress_logger.debug("Solving DR polishing problem") - # NOTE: this objective evalaution may not be accurate, due - # to the current initialization scheme for the auxiliary - # variables. new initialization will be implemented in the - # near future. - polishing_obj = polishing_model.scenarios[0, 0].polishing_obj + # polishing objective should be consistent with value of sum + # of absolute values of polynomial DR terms provided + # auxiliary variables initialized correctly + polishing_obj = polishing_model.polishing_obj config.progress_logger.debug(f" Initial DR norm: {value(polishing_obj)}") # === Solve the polishing model @@ -533,14 +490,14 @@ def minimize_dr_vars(model_data, config): # purposes config.progress_logger.debug(" Done solving DR polishing problem") config.progress_logger.debug( - f" Solve time: {getattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR)} s" + f" Termination condition: {results.solver.termination_condition} " ) config.progress_logger.debug( - f" Termination status: {results.solver.termination_condition} " + f" Solve time: {getattr(results.solver, TIC_TOC_SOLVE_TIME_ATTR)} s" ) # === Process solution by termination condition - acceptable = {tc.globallyOptimal, tc.optimal, tc.locallyOptimal, tc.feasible} + acceptable = {tc.globallyOptimal, tc.optimal, tc.locallyOptimal} if results.solver.termination_condition not in acceptable: # continue with "unpolished" master model solution config.progress_logger.warning( @@ -556,6 +513,7 @@ def minimize_dr_vars(model_data, config): # update master model second-stage, state, and decision rule # variables to polishing model solution polishing_model.solutions.load_from(results) + for idx, blk in model_data.master_model.scenarios.items(): ssv_zip = zip( blk.util.second_stage_variables, @@ -579,29 +537,29 @@ def minimize_dr_vars(model_data, config): mvar.set_value(value(pvar), skip_validation=True) config.progress_logger.debug(f" Optimized DR norm: {value(polishing_obj)}") - config.progress_logger.debug(" Polished Master objective:") + config.progress_logger.debug(" Polished master objective:") - # print master solution + # print breakdown of objective value of polished master solution if config.objective_focus == ObjectiveType.worst_case: - worst_blk_idx = max( + eval_obj_blk_idx = max( model_data.master_model.scenarios.keys(), key=lambda idx: value( model_data.master_model.scenarios[idx].second_stage_objective ), ) else: - worst_blk_idx = (0, 0) + eval_obj_blk_idx = (0, 0) # debugging: summarize objective breakdown - worst_master_blk = model_data.master_model.scenarios[worst_blk_idx] + eval_obj_blk = model_data.master_model.scenarios[eval_obj_blk_idx] config.progress_logger.debug( - " First-stage objective: " f"{value(worst_master_blk.first_stage_objective)}" + " First-stage objective: " f"{value(eval_obj_blk.first_stage_objective)}" ) config.progress_logger.debug( - " Second-stage objective: " f"{value(worst_master_blk.second_stage_objective)}" + " Second-stage objective: " f"{value(eval_obj_blk.second_stage_objective)}" ) polished_master_obj = value( - worst_master_blk.first_stage_objective + worst_master_blk.second_stage_objective + eval_obj_blk.first_stage_objective + eval_obj_blk.second_stage_objective ) config.progress_logger.debug(f" Objective: {polished_master_obj}") @@ -651,37 +609,64 @@ def add_scenario_to_master(model_data, violations): return -def higher_order_decision_rule_efficiency(config, model_data): - # === Efficiencies for decision rules - # if iteration <= |q| then all d^n where n > 1 are fixed to 0 - # if iteration == 0, all d^n, n > 0 are fixed to 0 - # These efficiencies should be carried through as d* to polishing - nlp_model = model_data.master_model - if config.decision_rule_order != None and len(config.second_stage_variables) > 0: - # Ensure all are unfixed unless next conditions are met... - for dr_var in nlp_model.scenarios[0, 0].util.decision_rule_vars: - dr_var.unfix() - num_dr_vars = len( - nlp_model.scenarios[0, 0].util.decision_rule_vars[0] - ) # there is at least one dr var - num_uncertain_params = len(config.uncertain_params) - nlp_model.const_efficiency_applied = False - nlp_model.linear_efficiency_applied = False - if model_data.iteration == 0: - nlp_model.const_efficiency_applied = True - for dr_var in nlp_model.scenarios[0, 0].util.decision_rule_vars: - for i in range(1, num_dr_vars): - dr_var[i].fix(0) - elif ( - model_data.iteration <= num_uncertain_params - and config.decision_rule_order > 1 - ): - # Only applied in DR order > 1 case - for dr_var in nlp_model.scenarios[0, 0].util.decision_rule_vars: - for i in range(num_uncertain_params + 1, num_dr_vars): - nlp_model.linear_efficiency_applied = True - dr_var[i].fix(0) - return +def get_master_dr_degree(model_data, config): + """ + Determine DR polynomial degree to enforce based on + the iteration number. + + Currently, 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 + - min(2, config.decision_rule_order) otherwise. + + Parameters + ---------- + model_data : MasterProblemData + Master problem data. + config : ConfigDict + PyROS solver options. + + Returns + ------- + int + DR order, or polynomial degree, to enforce. + """ + if model_data.iteration == 0: + return 0 + elif model_data.iteration <= len(config.uncertain_params): + return min(1, config.decision_rule_order) + else: + return min(2, config.decision_rule_order) + + +def higher_order_decision_rule_efficiency(model_data, config): + """ + Enforce DR coefficient variable efficiencies for + master problem-like formulation. + + Parameters + ---------- + model_data : MasterProblemData + Master problem data. + config : ConfigDict + PyROS solver options. + + Note + ---- + The DR coefficient variable efficiencies consist of + setting the degree of the DR polynomial expressions + by fixing the appropriate variables to 0. The degree + to be set depends on the iteration number; + see ``get_master_dr_degree``. + """ + order_to_enforce = get_master_dr_degree(model_data, config) + enforce_dr_degree( + blk=model_data.master_model.scenarios[0, 0], + config=config, + degree=order_to_enforce, + ) def solver_call_master(model_data, config, solver, solve_data): @@ -717,7 +702,7 @@ def solver_call_master(model_data, config, solver, solve_data): else: solvers = [solver] + config.backup_local_solvers - higher_order_decision_rule_efficiency(config, model_data) + higher_order_decision_rule_efficiency(model_data=model_data, config=config) solve_mode = "global" if config.solve_master_globally else "local" config.progress_logger.debug("Solving master problem") @@ -812,17 +797,28 @@ def solver_call_master(model_data, config, solver, solve_data): ) # debugging: log breakdown of master objective + if config.objective_focus == ObjectiveType.worst_case: + eval_obj_blk_idx = max( + nlp_model.scenarios.keys(), + key=lambda idx: value( + nlp_model.scenarios[idx].second_stage_objective + ), + ) + else: + eval_obj_blk_idx = (0, 0) + + eval_obj_blk = nlp_model.scenarios[eval_obj_blk_idx] config.progress_logger.debug(" Optimized master objective breakdown:") config.progress_logger.debug( - f" First-stage objective: {master_soln.first_stage_objective}" + f" First-stage objective: {value(eval_obj_blk.first_stage_objective)}" ) config.progress_logger.debug( - f" Second-stage objective: {master_soln.second_stage_objective}" + f" Second-stage objective: {value(eval_obj_blk.second_stage_objective)}" ) master_obj = ( - master_soln.first_stage_objective + master_soln.second_stage_objective + eval_obj_blk.first_stage_objective + eval_obj_blk.second_stage_objective ) - config.progress_logger.debug(f" Objective: {master_obj}") + config.progress_logger.debug(f" Objective: {value(master_obj)}") config.progress_logger.debug( f" Termination condition: {results.solver.termination_condition}" ) diff --git a/pyomo/contrib/pyros/pyros.py b/pyomo/contrib/pyros/pyros.py index 5b37b114722..829184fc70c 100644 --- a/pyomo/contrib/pyros/pyros.py +++ b/pyomo/contrib/pyros/pyros.py @@ -49,7 +49,7 @@ from datetime import datetime -__version__ = "1.2.8" +__version__ = "1.2.9" default_pyros_solver_logger = setup_pyros_logger() @@ -990,13 +990,14 @@ def solve( # === Move bounds on control variables to explicit ineq constraints wm_util = model_data.working_model - # === Assuming all other Var objects in the model are state variables + # === Every non-fixed variable that is neither first-stage + # nor second-stage is taken to be a state variable fsv = ComponentSet(model_data.working_model.util.first_stage_variables) ssv = ComponentSet(model_data.working_model.util.second_stage_variables) sv = ComponentSet() model_data.working_model.util.state_vars = [] for v in model_data.working_model.component_data_objects(Var): - if v not in fsv and v not in ssv and v not in sv: + if not v.fixed and v not in fsv | ssv | sv: model_data.working_model.util.state_vars.append(v) sv.add(v) diff --git a/pyomo/contrib/pyros/pyros_algorithm_methods.py b/pyomo/contrib/pyros/pyros_algorithm_methods.py index af7a91d21a4..38f675e64a5 100644 --- a/pyomo/contrib/pyros/pyros_algorithm_methods.py +++ b/pyomo/contrib/pyros/pyros_algorithm_methods.py @@ -886,6 +886,16 @@ def ROSolver_iterative_solve(model_data, config): np.array([pt for pt in separation_data.points_added_to_master]) ) + # initialize second-stage and state variables + # for new master block to separation + # solution chosen by heuristic. consequently, + # equality constraints should all be satisfied (up to tolerances). + for var, val in separation_results.violating_separation_variable_values.items(): + master_var = master_data.master_model.scenarios[k + 1, 0].find_component( + var + ) + master_var.set_value(val) + k += 1 iter_log_record.log(config.progress_logger.info) diff --git a/pyomo/contrib/pyros/separation_problem_methods.py b/pyomo/contrib/pyros/separation_problem_methods.py index 61f347b418d..240291f5375 100644 --- a/pyomo/contrib/pyros/separation_problem_methods.py +++ b/pyomo/contrib/pyros/separation_problem_methods.py @@ -882,13 +882,33 @@ def initialize_separation(perf_con_to_maximize, model_data, config): discrete geometry (as there is no master model block corresponding to any of the remaining discrete scenarios against which we separate). + + This method assumes that the master model has only one block + per iteration. """ - # initialize to values from nominal block if nominal objective. - # else, initialize to values from latest block added to master - if config.objective_focus == ObjectiveType.nominal: - block_num = 0 - else: - block_num = model_data.iteration + + def eval_master_violation(block_idx): + """ + Evaluate violation of `perf_con` by variables of + specified master block. + """ + new_con_map = ( + model_data.separation_model.util.map_new_constraint_list_to_original_con + ) + in_new_cons = perf_con_to_maximize in new_con_map + if in_new_cons: + sep_con = new_con_map[perf_con_to_maximize] + else: + sep_con = perf_con_to_maximize + master_con = model_data.master_model.scenarios[block_idx, 0].find_component( + sep_con + ) + return value(master_con) + + # initialize from master block with max violation of the + # performance constraint of interest. This gives the best known + # feasible solution (for case of non-discrete uncertainty sets). + block_num = max(range(model_data.iteration + 1), key=eval_master_violation) master_blk = model_data.master_model.scenarios[block_num, 0] master_blks = list(model_data.master_model.scenarios.values()) diff --git a/pyomo/contrib/pyros/tests/test_grcs.py b/pyomo/contrib/pyros/tests/test_grcs.py index 46af1277ba5..05cbdb849f4 100644 --- a/pyomo/contrib/pyros/tests/test_grcs.py +++ b/pyomo/contrib/pyros/tests/test_grcs.py @@ -5,10 +5,15 @@ import pyomo.common.unittest as unittest from pyomo.common.log import LoggingIntercept -from pyomo.common.collections import ComponentSet +from pyomo.common.collections import ComponentSet, ComponentMap from pyomo.common.config import ConfigBlock, ConfigValue from pyomo.core.base.set_types import NonNegativeIntegers -from pyomo.core.expr import identify_variables, identify_mutable_parameters +from pyomo.core.expr import ( + identify_variables, + identify_mutable_parameters, + MonomialTermExpression, + SumExpression, +) from pyomo.contrib.pyros.util import ( selective_clone, add_decision_rule_variables, @@ -27,8 +32,21 @@ from pyomo.contrib.pyros.util import identify_objective_functions from pyomo.common.collections import Bunch import time +import math from pyomo.contrib.pyros.util import time_code -from pyomo.contrib.pyros.uncertainty_sets import * +from pyomo.contrib.pyros.uncertainty_sets import ( + UncertaintySet, + BoxSet, + CardinalitySet, + BudgetSet, + FactorModelSet, + PolyhedralSet, + EllipsoidalSet, + AxisAlignedEllipsoidalSet, + IntersectionSet, + DiscreteScenarioSet, + Geometry, +) from pyomo.contrib.pyros.master_problem_methods import ( add_scenario_to_master, initial_construct_master, @@ -48,6 +66,11 @@ Solution, ) from pyomo.environ import ( + Reals, + Set, + Block, + ConstraintList, + ConcreteModel, Constraint, Expression, Objective, @@ -60,8 +83,11 @@ sin, sqrt, value, + maximize, + minimize, ) import logging +from itertools import chain logger = logging.getLogger(__name__) @@ -237,183 +263,351 @@ def test_cloning_positive_case(self): class testAddDecisionRuleVars(unittest.TestCase): - ''' - Testing the method to add decision rule variables to a Pyomo model. This function should add decision rule - variables to the list of first_stage_variables in a model object. The number of decision rule variables added - depends on the number of control variables in the model and the number of uncertain parameters in the model. - ''' + """ + Test method for adding decision rule variables to working model. + The number of decision rule variables per control variable + should depend on: - @unittest.skipIf(not scipy_available, 'Scipy is not available.') - def test_add_decision_rule_vars_positive_case(self): - ''' - Testing whether the correct number of decision rule variables is created in each DR type case - ''' + - the number of uncertain parameters in the model + - the decision rule order specified by the user. + """ + + def make_simple_test_model(self): + """ + Make simple test model for DR variable + declaration testing. + """ m = ConcreteModel() - m.p1 = Param(initialize=0, mutable=True) - m.p2 = Param(initialize=0, mutable=True) - m.z1 = Var(initialize=0) - m.z2 = Var(initialize=0) - m.working_model = ConcreteModel() - m.working_model.util = Block() + # uncertain parameters + m.p = Param(range(3), initialize=0, mutable=True) - m.working_model.util.second_stage_variables = [m.z1, m.z2] - m.working_model.util.uncertain_params = [m.p1, m.p2] - m.working_model.util.first_stage_variables = [] + # second-stage variables + m.z = Var([0, 1], initialize=0) - m.working_model.util.first_stage_variables = [] - config = Block() + # util block + m.util = Block() + m.util.first_stage_variables = [] + m.util.second_stage_variables = list(m.z.values()) + m.util.uncertain_params = list(m.p.values()) + + return m + + @unittest.skipIf(not scipy_available, 'Scipy is not available.') + def test_correct_num_dr_vars_static(self): + """ + Test DR variable setup routines declare the correct + number of DR coefficient variables, static DR case. + """ + model_data = ROSolveResults() + model_data.working_model = m = self.make_simple_test_model() + + config = Bunch() config.decision_rule_order = 0 - add_decision_rule_variables(model_data=m, config=config) + add_decision_rule_variables(model_data=model_data, config=config) + + for indexed_dr_var in m.util.decision_rule_vars: + self.assertEqual( + len(indexed_dr_var), + 1, + msg=( + "Number of decision rule coefficient variables " + f"in indexed Var object {indexed_dr_var.name!r}" + "does not match correct value." + ), + ) self.assertEqual( - len(m.working_model.util.first_stage_variables), - len(m.working_model.util.second_stage_variables), - msg="For static approximation decision rule the number of decision rule variables" - "added to the list of design variables should equal the number of control variables.", + len(ComponentSet(m.util.decision_rule_vars)), + len(m.util.second_stage_variables), + msg=( + "Number of unique indexed DR variable components should equal " + "number of second-stage variables." + ), ) - m.working_model.util.first_stage_variables = [] - - m.working_model.del_component(m.working_model.decision_rule_var_0) - m.working_model.del_component(m.working_model.decision_rule_var_1) + @unittest.skipIf(not scipy_available, 'Scipy is not available.') + def test_correct_num_dr_vars_affine(self): + """ + Test DR variable setup routines declare the correct + number of DR coefficient variables, affine DR case. + """ + model_data = ROSolveResults() + model_data.working_model = m = self.make_simple_test_model() + config = Bunch() config.decision_rule_order = 1 - add_decision_rule_variables(m, config=config) + add_decision_rule_variables(model_data=model_data, config=config) + + for indexed_dr_var in m.util.decision_rule_vars: + self.assertEqual( + len(indexed_dr_var), + 1 + len(m.util.uncertain_params), + msg=( + "Number of decision rule coefficient variables " + f"in indexed Var object {indexed_dr_var.name!r}" + "does not match correct value." + ), + ) self.assertEqual( - len(m.working_model.util.first_stage_variables), - len(m.working_model.util.second_stage_variables) - * (1 + len(m.working_model.util.uncertain_params)), - msg="For affine decision rule the number of decision rule variables add to the " - "list of design variables should equal the number of control variables" - "multiplied by the number of uncertain parameters plus 1.", + len(ComponentSet(m.util.decision_rule_vars)), + len(m.util.second_stage_variables), + msg=( + "Number of unique indexed DR variable components should equal " + "number of second-stage variables." + ), ) - m.working_model.util.first_stage_variables = [] - - m.working_model.del_component(m.working_model.decision_rule_var_0) - m.working_model.del_component(m.working_model.decision_rule_var_1) - m.working_model.del_component(m.working_model.decision_rule_var_0_index) - m.working_model.del_component(m.working_model.decision_rule_var_1_index) + @unittest.skipIf(not scipy_available, 'Scipy is not available.') + def test_correct_num_dr_vars_quadratic(self): + """ + Test DR variable setup routines declare the correct + number of DR coefficient variables, quadratic DR case. + """ + model_data = ROSolveResults() + model_data.working_model = m = self.make_simple_test_model() + config = Bunch() config.decision_rule_order = 2 - add_decision_rule_variables(m, config=config) + add_decision_rule_variables(model_data=model_data, config=config) + + num_params = len(m.util.uncertain_params) + correct_num_dr_vars = ( + 1 # static term + + num_params # affine terms + + sp.special.comb(num_params, 2, repetition=True, exact=True) + # quadratic terms + ) + for indexed_dr_var in m.util.decision_rule_vars: + self.assertEqual( + len(indexed_dr_var), + correct_num_dr_vars, + msg=( + "Number of decision rule coefficient variables " + f"in indexed Var object {indexed_dr_var.name!r}" + "does not match correct value." + ), + ) self.assertEqual( - len(m.working_model.util.first_stage_variables), - len(m.working_model.util.second_stage_variables) - * int( - 2 * len(m.working_model.util.uncertain_params) - + sp.special.comb(N=len(m.working_model.util.uncertain_params), k=2) - + 1 + len(ComponentSet(m.util.decision_rule_vars)), + len(m.util.second_stage_variables), + msg=( + "Number of unique indexed DR variable components should equal " + "number of second-stage variables." ), - msg="For quadratic decision rule the number of decision rule variables add to the " - "list of design variables should equal the number of control variables" - "multiplied by 2 time the number of uncertain parameters plus all 2-combinations" - "of uncertain parameters plus 1.", ) class testAddDecisionRuleConstraints(unittest.TestCase): - ''' - Testing the addition of decision rule constraints functionally relating second-stage (control) variables to - uncertain parameters and decision rule variables. This method should add constraints to the model object equal - to the number of control variables. These constraints should reference the uncertain parameters and unique - decision rule variables per control variable. - ''' + """ + Test method for adding decision rule equality constraints + to the working model. There should be as many decision + rule equality constraints as there are second-stage + variables, and each constraint should relate a second-stage + variable to the uncertain parameters and corresponding + decision rule variables. + """ - def test_correct_number_of_decision_rule_constraints(self): - ''' - Number of decision rule constraints added to the model should equal number of control variables in - list "second_stage_variables". - ''' + def make_simple_test_model(self): + """ + Make simple model for DR constraint testing. + """ m = ConcreteModel() - m.p1 = Param(initialize=0, mutable=True) - m.p2 = Param(initialize=0, mutable=True) - m.z1 = Var(initialize=0) - m.z2 = Var(initialize=0) - m.working_model = ConcreteModel() - m.working_model.util = Block() + # uncertain parameters + m.p = Param(range(3), initialize=0, mutable=True) - # === Decision rule vars have been added - m.working_model.decision_rule_var_0 = Var(initialize=0) - m.working_model.decision_rule_var_1 = Var(initialize=0) + # second-stage variables + m.z = Var([0, 1], initialize=0) - m.working_model.util.second_stage_variables = [m.z1, m.z2] - m.working_model.util.uncertain_params = [m.p1, m.p2] + # util block + m.util = Block() + m.util.first_stage_variables = [] + m.util.second_stage_variables = list(m.z.values()) + m.util.uncertain_params = list(m.p.values()) - decision_rule_cons = [] - config = Block() - config.decision_rule_order = 0 + return m - add_decision_rule_constraints(model_data=m, config=config) + @unittest.skipIf(not scipy_available, 'Scipy is not available.') + def test_num_dr_eqns_added_correct(self): + """ + Check that number of DR equality constraints added + by constraint declaration routines matches the number + of second-stage variables in the model. + """ + model_data = ROSolveResults() + model_data.working_model = m = self.make_simple_test_model() + + # === Decision rule vars have been added + m.decision_rule_var_0 = Var([0], initialize=0) + m.decision_rule_var_1 = Var([0], initialize=0) + m.util.decision_rule_vars = [m.decision_rule_var_0, m.decision_rule_var_1] + + # set up simple config-like object + config = Bunch() + config.decision_rule_order = 0 - for c in m.working_model.component_data_objects(Constraint, descend_into=True): - if "decision_rule_eqn_" in c.name: - decision_rule_cons.append(c) - m.working_model.del_component(c) + add_decision_rule_constraints(model_data=model_data, config=config) self.assertEqual( - len(decision_rule_cons), - len(m.working_model.util.second_stage_variables), + len(m.util.decision_rule_eqns), + len(m.util.second_stage_variables), msg="The number of decision rule constraints added to model should equal" "the number of control variables in the model.", ) - decision_rule_cons = [] - config.decision_rule_order = 1 + @unittest.skipIf(not scipy_available, 'Scipy is not available.') + def test_dr_eqns_form_correct(self): + """ + Check that form of decision rule equality constraints + is as expected. - # === Decision rule vars have been added - m.working_model.del_component(m.working_model.decision_rule_var_0) - m.working_model.del_component(m.working_model.decision_rule_var_1) + Decision rule equations should be of the standard form: + (sum of DR monomial terms) - (second-stage variable) == 0 + where each monomial term should be of form: + (product of uncertain parameters) * (decision rule variable) - m.working_model.decision_rule_var_0 = Var([0, 1, 2], initialize=0) - m.working_model.decision_rule_var_1 = Var([0, 1, 2], initialize=0) + This test checks that the equality constraints are of this + standard form. + """ + # set up simple model data like object + model_data = ROSolveResults() + model_data.working_model = m = self.make_simple_test_model() - add_decision_rule_constraints(model_data=m, config=config) + # set up simple config-like object + config = Bunch() + config.decision_rule_order = 2 - for c in m.working_model.component_data_objects(Constraint, descend_into=True): - if "decision_rule_eqn_" in c.name: - decision_rule_cons.append(c) - m.working_model.del_component(c) + # add DR variables and constraints + add_decision_rule_variables(model_data, config) + add_decision_rule_constraints(model_data, config) + + # DR polynomial terms and order in which they should + # appear depends on number of uncertain parameters + # and order in which the parameters are listed. + # so uncertain parameters participating in each term + # of the monomial is known, and listed out here. + dr_monomial_param_combos = [ + (1,), + (m.p[0],), + (m.p[1],), + (m.p[2],), + (m.p[0], m.p[0]), + (m.p[0], m.p[1]), + (m.p[0], m.p[2]), + (m.p[1], m.p[1]), + (m.p[1], m.p[2]), + (m.p[2], m.p[2]), + ] - self.assertEqual( - len(decision_rule_cons), - len(m.working_model.util.second_stage_variables), - msg="The number of decision rule constraints added to model should equal" - "the number of control variables in the model.", + dr_zip = zip( + m.util.second_stage_variables, + m.util.decision_rule_vars, + m.util.decision_rule_eqns, ) + for ss_var, indexed_dr_var, dr_eq in dr_zip: + dr_eq_terms = dr_eq.body.args - decision_rule_cons = [] - config.decision_rule_order = 2 + # check constraint body is sum expression + self.assertTrue( + isinstance(dr_eq.body, SumExpression), + msg=( + f"Body of DR constraint {dr_eq.name!r} is not of type " + f"{SumExpression.__name__}." + ), + ) - # === Decision rule vars have been added - m.working_model.del_component(m.working_model.decision_rule_var_0) - m.working_model.del_component(m.working_model.decision_rule_var_1) - m.working_model.del_component(m.working_model.decision_rule_var_0_index) - m.working_model.del_component(m.working_model.decision_rule_var_1_index) + # ensure DR equation has correct number of (additive) terms + self.assertEqual( + len(dr_eq_terms), + len(dr_monomial_param_combos) + 1, + msg=( + "Number of additive terms in the DR expression of " + f"DR constraint with name {dr_eq.name!r} does not match " + "expected value." + ), + ) - m.working_model.decision_rule_var_0 = Var([0, 1, 2, 3, 4, 5], initialize=0) - m.working_model.decision_rule_var_1 = Var([0, 1, 2, 3, 4, 5], initialize=0) + # check last term is negative of second-stage variable + second_stage_var_term = dr_eq_terms[-1] + last_term_is_neg_ss_var = ( + isinstance(second_stage_var_term, MonomialTermExpression) + and (second_stage_var_term.args[0] == -1) + and (second_stage_var_term.args[1] is ss_var) + and len(second_stage_var_term.args) == 2 + ) + self.assertTrue( + last_term_is_neg_ss_var, + msg=( + "Last argument of last term in second-stage variable" + f"term of DR constraint with name {dr_eq.name!r} " + "is not the negative corresponding second-stage variable " + f"{ss_var.name!r}" + ), + ) - add_decision_rule_constraints(model_data=m, config=config) + # now we check the other terms. + # these should comprise the DR polynomial expression + dr_polynomial_terms = dr_eq_terms[:-1] + dr_polynomial_zip = zip( + dr_polynomial_terms, indexed_dr_var.values(), dr_monomial_param_combos + ) + for idx, (term, dr_var, param_combo) in enumerate(dr_polynomial_zip): + # term should be a monomial expression of form + # (uncertain parameter product) * (decision rule variable) + # so length of expression object should be 2 + self.assertEqual( + len(term.args), + 2, + msg=( + f"Length of `args` attribute of term {str(term)} " + f"of DR equation {dr_eq.name!r} is not as expected. " + f"Args: {term.args}" + ), + ) - for c in m.working_model.component_data_objects(Constraint, descend_into=True): - if "decision_rule_eqn_" in c.name: - decision_rule_cons.append(c) - m.working_model.del_component(c) + # check that uncertain parameters participating in + # the monomial are as expected + param_product_multiplicand = term.args[0] + if idx == 0: + # static DR term + param_combo_found_in_term = (param_product_multiplicand,) + param_names = (str(param) for param in param_combo) + elif len(param_combo) == 1: + # affine DR terms + param_combo_found_in_term = (param_product_multiplicand,) + param_names = (param.name for param in param_combo) + else: + # higher-order DR terms + param_combo_found_in_term = param_product_multiplicand.args + param_names = (param.name for param in param_combo) + + self.assertEqual( + param_combo_found_in_term, + param_combo, + msg=( + f"All but last multiplicand of DR monomial {str(term)} " + f"is not the uncertain parameter tuple " + f"({', '.join(param_names)})." + ), + ) - self.assertEqual( - len(decision_rule_cons), - len(m.working_model.util.second_stage_variables), - msg="The number of decision rule constraints added to model should equal" - "the number of control variables in the model.", - ) + # check that DR variable participating in the monomial + # is as expected + dr_var_multiplicand = term.args[1] + self.assertIs( + dr_var_multiplicand, + dr_var, + msg=( + f"Last multiplicand of DR monomial {str(term)} " + f"is not the DR variable {dr_var.name!r}." + ), + ) class testModelIsValid(unittest.TestCase): @@ -3572,6 +3766,9 @@ def test_solve_master(self): master_data.master_model.scenarios[0, 0].second_stage_objective = Expression( expr=master_data.master_model.scenarios[0, 0].x ) + master_data.master_model.scenarios[ + 0, 0 + ].util.dr_var_to_exponent_map = ComponentMap() master_data.iteration = 0 master_data.timing = TimingData() diff --git a/pyomo/contrib/pyros/util.py b/pyomo/contrib/pyros/util.py index 19f178c70f6..6aa6ed61936 100644 --- a/pyomo/contrib/pyros/util.py +++ b/pyomo/contrib/pyros/util.py @@ -3,7 +3,7 @@ ''' import copy from enum import Enum, auto -from pyomo.common.collections import ComponentSet +from pyomo.common.collections import ComponentSet, ComponentMap from pyomo.common.modeling import unique_component_name from pyomo.core.base import ( Constraint, @@ -17,11 +17,11 @@ Block, Param, ) +from pyomo.core.util import prod from pyomo.core.base.var import IndexedVar from pyomo.core.base.set_types import Reals from pyomo.opt import TerminationCondition as tc from pyomo.core.expr import value -import pyomo.core.expr as EXPR from pyomo.core.expr.numeric_expr import NPV_MaxExpression, NPV_MinExpression from pyomo.repn.standard_repn import generate_standard_repn from pyomo.core.expr.visitor import ( @@ -39,7 +39,6 @@ import timeit from contextlib import contextmanager import logging -from pprint import pprint import math from pyomo.common.timing import HierarchicalTimer from pyomo.common.log import Preformatted @@ -1275,220 +1274,149 @@ def selective_clone(block, first_stage_vars): def add_decision_rule_variables(model_data, config): - ''' - Function to add decision rule (DR) variables to the working model. DR variables become first-stage design - variables which do not get copied at each iteration. Currently support static_approx (no DR), affine DR, - and quadratic DR. - :param model_data: the data container for the working model - :param config: the config block - :return: - ''' + """ + Add variables for polynomial decision rules to the working + model. + + Parameters + ---------- + model_data : ROSolveResults + Model data. + config : config_dict + PyROS solver options. + + Note + ---- + Decision rule variables are considered first-stage decision + variables which do not get copied at each iteration. + PyROS currently supports static (zeroth order), + affine (first-order), and quadratic DR. + """ second_stage_variables = model_data.working_model.util.second_stage_variables first_stage_variables = model_data.working_model.util.first_stage_variables - uncertain_params = model_data.working_model.util.uncertain_params decision_rule_vars = [] + + # since DR expression is a general polynomial in the uncertain + # parameters, the exact number of DR variables per second-stage + # variable depends on DR order and uncertainty set dimension degree = config.decision_rule_order - bounds = (None, None) - if degree == 0: - for i in range(len(second_stage_variables)): - model_data.working_model.add_component( - "decision_rule_var_" + str(i), - Var( - initialize=value(second_stage_variables[i], exception=False), - bounds=bounds, - domain=Reals, - ), - ) - first_stage_variables.extend( - getattr( - model_data.working_model, "decision_rule_var_" + str(i) - ).values() - ) - decision_rule_vars.append( - getattr(model_data.working_model, "decision_rule_var_" + str(i)) - ) - elif degree == 1: - for i in range(len(second_stage_variables)): - index_set = list(range(len(uncertain_params) + 1)) - model_data.working_model.add_component( - "decision_rule_var_" + str(i), - Var(index_set, initialize=0, bounds=bounds, domain=Reals), - ) - # === For affine drs, the [0]th constant term is initialized to the control variable values, all other terms are initialized to 0 - getattr(model_data.working_model, "decision_rule_var_" + str(i))[ - 0 - ].set_value( - value(second_stage_variables[i], exception=False), skip_validation=True - ) - first_stage_variables.extend( - list( - getattr( - model_data.working_model, "decision_rule_var_" + str(i) - ).values() - ) - ) - decision_rule_vars.append( - getattr(model_data.working_model, "decision_rule_var_" + str(i)) - ) - elif degree == 2 or degree == 3 or degree == 4: - for i in range(len(second_stage_variables)): - num_vars = int(sp.special.comb(N=len(uncertain_params) + degree, k=degree)) - dict_init = {} - for r in range(num_vars): - if r == 0: - dict_init.update( - {r: value(second_stage_variables[i], exception=False)} - ) - else: - dict_init.update({r: 0}) - model_data.working_model.add_component( - "decision_rule_var_" + str(i), - Var( - list(range(num_vars)), - initialize=dict_init, - bounds=bounds, - domain=Reals, - ), - ) - first_stage_variables.extend( - list( - getattr( - model_data.working_model, "decision_rule_var_" + str(i) - ).values() - ) - ) - decision_rule_vars.append( - getattr(model_data.working_model, "decision_rule_var_" + str(i)) - ) - else: - raise ValueError( - "Decision rule order " - + str(config.decision_rule_order) - + " is not yet supported. PyROS supports polynomials of degree 0 (static approximation), 1, 2." + num_uncertain_params = len(model_data.working_model.util.uncertain_params) + num_dr_vars = sp.special.comb( + N=num_uncertain_params + degree, k=degree, exact=True, repetition=False + ) + + for idx, ss_var in enumerate(second_stage_variables): + # declare DR coefficients for current second-stage variable + indexed_dr_var = Var( + range(num_dr_vars), initialize=0, bounds=(None, None), domain=Reals + ) + model_data.working_model.add_component( + f"decision_rule_var_{idx}", indexed_dr_var ) - model_data.working_model.util.decision_rule_vars = decision_rule_vars + # index 0 entry of the IndexedVar is the static + # DR term. initialize to user-provided value of + # the corresponding second-stage variable. + # all other entries remain initialized to 0. + indexed_dr_var[0].set_value(value(ss_var, exception=False)) -def partition_powers(n, v): - """Partition a total degree n across v variables + # update attributes + first_stage_variables.extend(indexed_dr_var.values()) + decision_rule_vars.append(indexed_dr_var) - This is an implementation of the "stars and bars" algorithm from - combinatorial mathematics. + model_data.working_model.util.decision_rule_vars = decision_rule_vars - This partitions a "total integer degree" of n across v variables - such that each variable gets an integer degree >= 0. You can think - of this as dividing a set of n+v things into v groupings, with the - power for each v_i being 1 less than the number of things in the - i'th group (because the v is part of the group). It is therefore - sufficient to just get the v-1 starting points chosen from a list of - indices n+v long (the first starting point is fixed to be 0). +def add_decision_rule_constraints(model_data, config): """ - for starts in it.combinations(range(1, n + v), v - 1): - # add the initial starting point to the beginning and the total - # number of objects (degree counters and variables) to the end - # of the list. The degree for each variable is 1 less than the - # difference of sequential starting points (to account for the - # variable itself) - starts = (0,) + starts + (n + v,) - yield [starts[i + 1] - starts[i] - 1 for i in range(v)] - - -def sort_partitioned_powers(powers_list): - powers_list = sorted(powers_list, reverse=True) - powers_list = sorted(powers_list, key=lambda elem: max(elem)) - return powers_list + Add decision rule equality constraints to the working model. - -def add_decision_rule_constraints(model_data, config): - ''' - Function to add the defining Constraint relationships for the decision rules to the working model. - :param model_data: model data container object - :param config: the config object - :return: - ''' + Parameters + ---------- + model_data : ROSolveResults + Model data. + config : ConfigDict + PyROS solver options. + """ second_stage_variables = model_data.working_model.util.second_stage_variables uncertain_params = model_data.working_model.util.uncertain_params decision_rule_eqns = [] + decision_rule_vars_list = model_data.working_model.util.decision_rule_vars degree = config.decision_rule_order - if degree == 0: - for i in range(len(second_stage_variables)): - model_data.working_model.add_component( - "decision_rule_eqn_" + str(i), - Constraint( - expr=getattr( - model_data.working_model, "decision_rule_var_" + str(i) - ) - == second_stage_variables[i] - ), - ) - decision_rule_eqns.append( - getattr(model_data.working_model, "decision_rule_eqn_" + str(i)) - ) - elif degree == 1: - for i in range(len(second_stage_variables)): - expr = 0 - for j in range( - len(getattr(model_data.working_model, "decision_rule_var_" + str(i))) - ): - if j == 0: - expr += getattr( - model_data.working_model, "decision_rule_var_" + str(i) - )[j] - else: - expr += ( - getattr( - model_data.working_model, "decision_rule_var_" + str(i) - )[j] - * uncertain_params[j - 1] - ) - model_data.working_model.add_component( - "decision_rule_eqn_" + str(i), - Constraint(expr=expr == second_stage_variables[i]), - ) - decision_rule_eqns.append( - getattr(model_data.working_model, "decision_rule_eqn_" + str(i)) - ) - elif degree >= 2: - # Using bars and stars groupings of variable powers, construct x1^a * .... * xn^b terms for all c <= a+...+b = degree - all_powers = [] - for n in range(1, degree + 1): - all_powers.append( - sort_partitioned_powers( - list(partition_powers(n, len(uncertain_params))) - ) - ) - for i in range(len(second_stage_variables)): - Z = list( - z - for z in getattr( - model_data.working_model, "decision_rule_var_" + str(i) - ).values() - ) - e = Z.pop(0) - for degree_param_powers in all_powers: - for param_powers in degree_param_powers: - product = 1 - for idx, power in enumerate(param_powers): - if power == 0: - pass - else: - product = product * uncertain_params[idx] ** power - e += Z.pop(0) * product - model_data.working_model.add_component( - "decision_rule_eqn_" + str(i), - Constraint(expr=e == second_stage_variables[i]), - ) - decision_rule_eqns.append( - getattr(model_data.working_model, "decision_rule_eqn_" + str(i)) + + # keeping track of degree of monomial in which each + # DR coefficient participates will be useful for later + dr_var_to_exponent_map = ComponentMap() + + # set up uncertain parameter combinations for + # construction of the monomials of the DR expressions + monomial_param_combos = [] + for power in range(degree + 1): + power_combos = it.combinations_with_replacement(uncertain_params, power) + monomial_param_combos.extend(power_combos) + + # now construct DR equations and declare them on the working model + second_stage_dr_var_zip = zip(second_stage_variables, decision_rule_vars_list) + for idx, (ss_var, indexed_dr_var) in enumerate(second_stage_dr_var_zip): + # for each DR equation, the number of coefficients should match + # the number of monomial terms exactly + if len(monomial_param_combos) != len(indexed_dr_var.index_set()): + raise ValueError( + f"Mismatch between number of DR coefficient variables " + f"and number of DR monomials for DR equation index {idx}, " + f"corresponding to second-stage variable {ss_var.name!r}. " + f"({len(indexed_dr_var.index_set())}!= {len(monomial_param_combos)})" ) - if len(Z) != 0: - raise RuntimeError( - "Construction of the decision rule functions did not work correctly! " - "Did not use all coefficient terms." - ) + + # construct the DR polynomial + dr_expression = 0 + for dr_var, param_combo in zip(indexed_dr_var.values(), monomial_param_combos): + dr_expression += dr_var * prod(param_combo) + + # map decision rule var to degree (exponent) of the + # associated monomial with respect to the uncertain params + dr_var_to_exponent_map[dr_var] = len(param_combo) + + # declare constraint on model + dr_eqn = Constraint(expr=dr_expression - ss_var == 0) + model_data.working_model.add_component(f"decision_rule_eqn_{idx}", dr_eqn) + + # append to list of DR equality constraints + decision_rule_eqns.append(dr_eqn) + + # finally, add attributes to util block model_data.working_model.util.decision_rule_eqns = decision_rule_eqns + model_data.working_model.util.dr_var_to_exponent_map = dr_var_to_exponent_map + + +def enforce_dr_degree(blk, config, degree): + """ + Make decision rule polynomials of a given degree + by fixing value of the appropriate subset of the decision + rule coefficients to 0. + + Parameters + ---------- + blk : ScalarBlock + Working model, or master problem block. + config : ConfigDict + PyROS solver options. + degree : int + Degree of the DR polynomials that is to be enforced. + """ + second_stage_vars = blk.util.second_stage_variables + indexed_dr_vars = blk.util.decision_rule_vars + dr_var_to_exponent_map = blk.util.dr_var_to_exponent_map + + for ss_var, indexed_dr_var in zip(second_stage_vars, indexed_dr_vars): + for dr_var in indexed_dr_var.values(): + dr_var_degree = dr_var_to_exponent_map[dr_var] + + if dr_var_degree > degree: + dr_var.fix(0) + else: + dr_var.unfix() def identify_objective_functions(model, objective):