# ___________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC # Under the terms of Contract DE-NA0003525 with National Technology and # Engineering Solutions of Sandia, LLC, the U.S. Government retains certain # rights in this software. # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ from pyomo.core import ConcreteModel, Var, Set, Objective, Constraint, ConstraintList, NonNegativeReals, NonPositiveReals, Expression, ComponentMap, maximize, minimize from pyomo.opt import (ProblemFormat, SolverFactory, SolverManagerFactory, SolverStatus, TerminationCondition, SolutionStatus) from pyomo.pysp.phutils import (isVariableNameIndexed, extractVariableNameAndIndex, extractComponentIndices) # # a routine to create the extensive form, given an input scenario tree and instances. # IMPT: unlike scenario instances, the extensive form instance is *not* self-contained. # in particular, it has binding constraints that cross the binding instance and # the scenario instances. it is up to the caller to keep track of which scenario # instances are associated with the extensive form. this might be something we # encapsulate at some later time. # # GAH: The following comment is still true, except that the scenario tree is no longer modified # # NOTE: if cvar terms are generated, then the input scenario tree is modified accordingly, # i.e., with the addition of the "eta" variable at the root node and the excess # variables at (for lack of a better place - they are per-scenario, but are not # blended) the second stage. # def create_ef_instance(scenario_tree, ef_instance_name = "MASTER", verbose_output = False, generate_weighted_cvar = False, cvar_weight = None, risk_alpha = None, cc_indicator_var_name=None, cc_alpha=0.0): # # create the new and empty binding instance. # # scenario tree must be "linked" with a set of instances # to used this function scenario_instances = {} for scenario in scenario_tree.scenarios: if scenario._instance is None: raise ValueError( "Cannot construct extensive form instance. " "The scenario tree does not appear to be linked " "to any Pyomo models. Missing model for scenario " "with name: %s" % (scenario.name)) scenario_instances[scenario.name] = scenario._instance binding_instance = ConcreteModel(name=ef_instance_name) root_node = scenario_tree.findRootNode() opt_sense = minimize \ if (scenario_tree._scenarios[0]._instance_objective.is_minimizing()) \ else maximize # # validate cvar options, if specified. # cvar_excess_vardatas = [] if generate_weighted_cvar: if (cvar_weight is None) or (cvar_weight < 0.0): raise RuntimeError("Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight)) if (risk_alpha is None) or (risk_alpha <= 0.0) or (risk_alpha >= 1.0): raise RuntimeError("CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha)) if verbose_output: print("Writing CVaR weighted objective") print("CVaR term weight="+str(cvar_weight)) print("CVaR alpha="+str(risk_alpha)) print("") # create the eta and excess variable on a per-scenario basis, # in addition to the constraint relating to the two. cvar_eta_variable_name = "CVAR_ETA_"+str(root_node._name) cvar_eta_variable = Var() binding_instance.add_component(cvar_eta_variable_name, cvar_eta_variable) excess_var_domain = NonNegativeReals if (opt_sense == minimize) else \ NonPositiveReals compute_excess_constraint = \ binding_instance.COMPUTE_SCENARIO_EXCESS = \ ConstraintList() for scenario in scenario_tree._scenarios: cvar_excess_variable_name = "CVAR_EXCESS_"+scenario._name cvar_excess_variable = Var(domain=excess_var_domain) binding_instance.add_component(cvar_excess_variable_name, cvar_excess_variable) compute_excess_expression = cvar_excess_variable compute_excess_expression -= scenario._instance_cost_expression compute_excess_expression += cvar_eta_variable if opt_sense == maximize: compute_excess_expression *= -1 compute_excess_constraint.add((0.0, compute_excess_expression, None)) cvar_excess_vardatas.append((cvar_excess_variable, scenario._probability)) # the individual scenario instances are sub-blocks of the binding instance. for scenario in scenario_tree._scenarios: scenario_instance = scenario_instances[scenario._name] binding_instance.add_component(str(scenario._name), scenario_instance) # Now deactivate the scenario instance Objective since we are creating # a new master objective scenario._instance_objective.deactivate() # walk the scenario tree - create variables representing the # common values for all scenarios associated with that node, along # with equality constraints to enforce non-anticipativity. also # create expected cost variables for each node, to be computed via # constraints/objectives defined in a subsequent pass. master # variables are created for all nodes but those in the last # stage. expected cost variables are, for no particularly good # reason other than easy coding, created for nodes in all stages. if verbose_output: print("Creating variables for master binding instance") _cmap = binding_instance.MASTER_CONSTRAINT_MAP = ComponentMap() for stage in scenario_tree._stages[:-1]: # skip the leaf stage for tree_node in stage._tree_nodes: # create the master blending variable and constraints for this node master_blend_variable_name = \ "MASTER_BLEND_VAR_"+str(tree_node._name) master_blend_constraint_name = \ "MASTER_BLEND_CONSTRAINT_"+str(tree_node._name) # don't create master variables for derived # stage variables as they will not be used in # the problem, and their values would likely # never be consistent with what is stored on the # scenario variables master_variable_index = Set(initialize=sorted(tree_node._standard_variable_ids), ordered=True, name=master_blend_variable_name+"_index") binding_instance.add_component(master_blend_variable_name+"_index", master_variable_index) master_variable = Var(master_variable_index, name=master_blend_variable_name) binding_instance.add_component(master_blend_variable_name, master_variable) master_constraint = ConstraintList(name=master_blend_constraint_name) binding_instance.add_component(master_blend_constraint_name, master_constraint) tree_node_variable_datas = tree_node._variable_datas for variable_id in sorted(tree_node._standard_variable_ids): master_vardata = master_variable[variable_id] vardatas = tree_node_variable_datas[variable_id] # Don't blend fixed variables if not tree_node.is_variable_fixed(variable_id): for scenario_vardata, scenario_probability in vardatas: _cmap[scenario_vardata] = master_constraint.add( (master_vardata-scenario_vardata,0.0)) if generate_weighted_cvar: cvar_cost_expression_name = "CVAR_COST_"+str(root_node._name) cvar_cost_expression = Expression(name=cvar_cost_expression_name) binding_instance.add_component(cvar_cost_expression_name, cvar_cost_expression) # create an expression to represent the expected cost at the root node binding_instance.EF_EXPECTED_COST = \ Expression(initialize=sum(scenario._probability * \ scenario._instance_cost_expression for scenario in scenario_tree._scenarios)) opt_expression = \ binding_instance.MASTER_OBJECTIVE_EXPRESSION = \ Expression(initialize=binding_instance.EF_EXPECTED_COST) if generate_weighted_cvar: cvar_cost_expression_name = "CVAR_COST_"+str(root_node._name) cvar_cost_expression = \ binding_instance.find_component(cvar_cost_expression_name) if cvar_weight == 0.0: # if the cvar weight is 0, then we're only # doing cvar - no mean. opt_expression.set_value(cvar_cost_expression) else: opt_expression.expr += cvar_weight * cvar_cost_expression binding_instance.MASTER = Objective(sense=opt_sense, expr=opt_expression) # CVaR requires the addition of a variable per scenario to # represent the cost excess, and a constraint to compute the cost # excess relative to eta. if generate_weighted_cvar: # add the constraint to compute the master CVaR variable value. iterate # over scenario instances to create the expected excess component first. cvar_cost_expression_name = "CVAR_COST_"+str(root_node._name) cvar_cost_expression = binding_instance.find_component(cvar_cost_expression_name) cvar_eta_variable_name = "CVAR_ETA_"+str(root_node._name) cvar_eta_variable = binding_instance.find_component(cvar_eta_variable_name) cost_expr = 1.0 for scenario_excess_vardata, scenario_probability in cvar_excess_vardatas: cost_expr += (scenario_probability * scenario_excess_vardata) cost_expr /= (1.0 - risk_alpha) cost_expr += cvar_eta_variable cvar_cost_expression.set_value(cost_expr) if cc_indicator_var_name is not None: if verbose_output is True: print("Creating chance constraint for indicator variable= "+cc_indicator_var_name) print( "with alpha= "+str(cc_alpha)) if not isVariableNameIndexed(cc_indicator_var_name): cc_expression = 0 #?????? for scenario in scenario_tree._scenarios: scenario_instance = scenario_instances[scenario._name] scenario_probability = scenario._probability cc_var = scenario_instance.find_component(cc_indicator_var_name) cc_expression += scenario_probability * cc_var def makeCCRule(expression): def CCrule(model): return(1.0 - cc_alpha, cc_expression, None) return CCrule cc_constraint_name = "cc_"+cc_indicator_var_name cc_constraint = Constraint(name=cc_constraint_name, rule=makeCCRule(cc_expression)) binding_instance.add_component(cc_constraint_name, cc_constraint) else: print("multiple cc not yet supported.") variable_name, index_template = extractVariableNameAndIndex(cc_indicator_var_name) # verify that the root variable exists and grab it. # NOTE: we are using whatever scenario happens to laying around... it might be better to use the reference variable = scenario_instance.find_component(variable_name) if variable is None: raise RuntimeError("Unknown variable="+variable_name+" referenced as the CC indicator variable.") # extract all "real", i.e., fully specified, indices matching the index template. match_indices = extractComponentIndices(variable, index_template) # there is a possibility that no indices match the input template. # if so, let the user know about it. if len(match_indices) == 0: raise RuntimeError("No indices match template="+str(index_template)+" for variable="+variable_name) # add the suffix to all variable values identified. for index in match_indices: variable_value = variable[index] cc_expression = 0 #?????? for scenario in scenario_tree._scenarios: scenario_instance = scenario_instances[scenario._name] scenario_probability = scenario._probability cc_var = scenario_instance.find_component(variable_name)[index] cc_expression += scenario_probability * cc_var def makeCCRule(expression): def CCrule(model): return(1.0 - cc_alpha, cc_expression, None) return CCrule indexasname = '' for c in str(index): if c not in ' ,': indexasname += c cc_constraint_name = "cc_"+variable_name+"_"+indexasname cc_constraint = Constraint(name=cc_constraint_name, rule=makeCCRule(cc_expression)) binding_instance.add_component(cc_constraint_name, cc_constraint) return binding_instance # # write the EF binding instance and all sub-instances. # def write_ef(binding_instance, output_filename, symbolic_solver_labels=False, output_fixed_variable_bounds=False): # determine the output file type, and invoke the appropriate # writer. pieces = output_filename.rsplit(".",1) if len(pieces) != 2: raise RuntimeError("Could not determine suffix from " "output filename="+output_filename) ef_output_file_suffix = pieces[1] io_options = {} if symbolic_solver_labels: io_options['symbolic_solver_labels'] = True if output_fixed_variable_bounds: io_options['output_fixed_variable_bounds'] = True # create the output file. if ef_output_file_suffix == "lp": _, smap_id = binding_instance.write( filename=output_filename, format=ProblemFormat.cpxlp, solver_capability=lambda x: True, io_options=io_options) elif ef_output_file_suffix == "nl": _, smap_id = binding_instance.write( filename=output_filename, format=ProblemFormat.nl, solver_capability=lambda x: True, io_options=io_options) elif ef_output_file_suffix == "mps": _, smap_id = binding_instance.write( filename=output_filename, format=ProblemFormat.mps, solver_capability=lambda x: True, io_options=io_options) else: raise RuntimeError("Unknown file suffix="+ef_output_file_suffix+ " specified when writing extensive form") return smap_id # # solve the EF binding instance and load the solution # def solve_ef(master_instance, options): with SolverFactory(options.solver_type) as ef_solver: with SolverManagerFactory(options.solver_manager_type) as ef_solver_manager: if ef_solver is None: raise ValueError("Failed to create solver manager of type=" +options.solver_type+" for use in extensive form solve") if len(options.ef_solver_options) > 0: print("Initializing ef solver with options="+str(options.ef_solver_options)) ef_solver.set_options("".join(options.ef_solver_options)) if options.ef_mipgap is not None: if (options.ef_mipgap < 0.0) or \ (options.ef_mipgap > 1.0): raise ValueError("Value of the mipgap parameter for the " "EF solve must be on the unit interval; " "value specified="+str(options.ef_mipgap)) else: ef_solver.mipgap = options.ef_mipgap solve_kwds = {} solve_kwds['load_solutions'] = False if options.keep_solver_files: solve_kwds['keepfiles'] = True if options.symbolic_solver_labels: solve_kwds['symbolic_solver_labels'] = True if options.output_solver_log: solve_kwds['tee'] = True if options.write_fixed_variables: solve_kwds['output_fixed_variable_bounds'] = True if options.verbose: print("Solving extensive form.") if (not options.disable_warmstarts) and \ (ef_solver.warm_start_capable()): ef_action_handle = ef_solver_manager.queue( master_instance, opt=ef_solver, warmstart=True, **solve_kwds) else: ef_action_handle = ef_solver_manager.queue( master_instance, opt=ef_solver, **solve_kwds) results = ef_solver_manager.wait_for(ef_action_handle) # check the return code - if this is anything but "have a solution", we need to bail. if (results.solver.status == SolverStatus.ok) and \ ((results.solver.termination_condition == TerminationCondition.optimal) or \ ((len(results.solution) > 0) and (results.solution(0).status == SolutionStatus.optimal))): master_instance.solutions.load_from(results) return results raise RuntimeError("Extensive form was infeasible!")