Changeset 2434
- Timestamp:
- Mar 12, 2010 3:43:03 AM (11 years ago)
- Location:
- coopr.pysp/stable/2.3
- Files:
-
- 11 edited
- 15 copied
Legend:
- Unmodified
- Added
- Removed
-
coopr.pysp/stable/2.3
- Property svnmerge-integrated changed
/coopr.pysp/trunk merged: 2393,2396,2398,2401-2406,2410-2414,2418,2423,2427-2428
- Property svnmerge-integrated changed
-
coopr.pysp/stable/2.3/coopr/pysp/__init__.py
r2317 r2434 21 21 from ef_writer_script import * 22 22 from phinit import * 23 from phsolvermanager import * 24 from phobjective import * 23 25 24 26 pyutilib.component.core.PluginGlobals.pop_env() 25 26 27 27 28 try: -
coopr.pysp/stable/2.3/coopr/pysp/ef.py
r2391 r2434 5 5 import traceback 6 6 import copy 7 import gc 7 8 8 9 from coopr.pysp.scenariotree import * … … 15 16 from coopr.pyomo.base.var import _VarValue, _VarBase 16 17 18 # 17 19 # brain-dead utility for determing if there is a binary to write in the 18 20 # composite model - need to know this, because CPLEX doesn't like empty 19 21 # binary blocks in the LP file. 22 # 23 20 24 def binaries_present(master_model, scenario_instances): 21 25 … … 34 38 return False 35 39 40 # 36 41 # brain-dead utility for determing if there is a binary to write in the 37 42 # composite model - need to know this, because CPLEX doesn't like empty 38 43 # integer blocks in the LP file. 44 # 45 39 46 def integers_present(master_model, scenario_instances): 40 47 … … 53 60 return False 54 61 62 # 63 # a routine to create the extensive form, given an input scenario tree and instances. 64 # IMPT: unlike scenario instances, the extensive form instance is *not* self-contained. 65 # in particular, it has binding constraints that cross the binding instance and 66 # the scenario instances. it is up to the caller to keep track of which scenario 67 # instances are associated with the extensive form. this might be something we 68 # encapsulate at some later time. 69 # 70 71 def create_ef_instance(scenario_tree, scenario_instances, 72 verbose_output = False, 73 generate_weighted_cvar = False, cvar_weight = None, risk_alpha = None): 74 75 # 76 # validate cvar options, if specified. 77 # 78 if generate_weighted_cvar is True: 79 if (cvar_weight is None) or (cvar_weight <= 0.0): 80 raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight) 81 if (risk_alpha is None) or (risk_alpha <= 0.0) or (risk_alpha >= 1.0): 82 raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha) 83 84 if verbose_output is True: 85 print "Writing CVaR weighted objective" 86 print "CVaR term weight="+str(cvar_weight) 87 print "CVaR alpha="+str(risk_alpha) 88 print "" 89 90 # 91 # create the new and empty binding instance. 92 # 93 94 binding_instance = Model() 95 binding_instance.name = "MASTER" 96 97 # walk the scenario tree - create variables representing the common values for all scenarios 98 # associated with that node. the constraints will be created later. also create expected-cost 99 # variables for each node, to be computed via constraints/objectives defined in a subsequent pass. 100 # master variables are created for all nodes but those in the last stage. expected cost variables 101 # are, for no particularly good reason other than easy coding, created for nodes in all stages. 102 if verbose_output is True: 103 print "Creating variables for master binding instance" 104 105 for stage in scenario_tree._stages: 106 107 # first loop is to create master (blended) variables across all stages but the last. 108 for (stage_variable, index_template, stage_variable_indices) in stage._variables: 109 110 if verbose_output is True: 111 print "Creating master variable and blending constraints for decision variable="+stage_variable.name+", indices="+str(index_template) 112 113 for tree_node in stage._tree_nodes: 114 115 if stage != scenario_tree._stages[-1]: 116 117 master_variable_name = tree_node._name + "_" + stage_variable.name 118 119 # because there may be a single stage variable and multiple indices, check 120 # for the existence of the variable at this node - if you don't, you'll 121 # inadvertently over-write what was there previously! 122 master_variable = None 123 try: 124 master_variable = getattr(binding_instance, master_variable_name) 125 except: 126 new_master_variable_index = stage_variable._index 127 new_master_variable = None 128 if (len(new_master_variable_index) is 1) and (None in new_master_variable_index): 129 new_master_variable = Var(name=stage_variable.name) 130 else: 131 new_master_variable = Var(new_master_variable_index, name=stage_variable.name) 132 new_master_variable.construct() 133 new_master_variable._model = binding_instance 134 setattr(binding_instance, master_variable_name, new_master_variable) 135 136 master_variable = new_master_variable 137 138 for index in stage_variable_indices: 139 140 is_used = True # until proven otherwise 141 for scenario in tree_node._scenarios: 142 instance = scenario_instances[scenario._name] 143 if getattr(instance,stage_variable.name)[index].status == VarStatus.unused: 144 is_used = False 145 146 is_fixed = False # until proven otherwise 147 for scenario in tree_node._scenarios: 148 instance = scenario_instances[scenario._name] 149 if getattr(instance,stage_variable.name)[index].fixed is True: 150 is_fixed = True 151 152 if (is_used is True) and (is_fixed is False): 153 154 # the following is necessary, specifically to get the name - deepcopy won't reset these attributes. 155 # and because presolve/simplification is name-based, the names *have* to be different. 156 master_variable[index].var = master_variable 157 master_variable[index].name = tree_node._name + "_" + master_variable[index].name 158 159 for scenario in tree_node._scenarios: 160 161 scenario_instance = scenario_instances[scenario._name] 162 scenario_variable = getattr(scenario_instance, stage_variable.name) 163 new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index) 164 new_constraint = Constraint(name=new_constraint_name) 165 new_expr = master_variable[index] - scenario_variable[index] 166 new_constraint.add(None, (0.0, new_expr, 0.0)) 167 new_constraint._model = binding_instance 168 setattr(binding_instance, new_constraint_name, new_constraint) 169 170 # the second loop is for creating the stage cost variable in each tree node. 171 for tree_node in stage._tree_nodes: 172 173 # create a variable to represent the expected cost at this node - 174 # the constraint to compute this comes later. 175 expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name 176 expected_cost_variable = Var(name=expected_cost_variable_name) 177 expected_cost_variable._model = binding_instance 178 setattr(binding_instance, expected_cost_variable_name, expected_cost_variable) 179 180 # if we're generating the weighted CVaR objective term, create the 181 # corresponding variable and the master CVaR eta variable. 182 if generate_weighted_cvar is True: 183 184 root_node = scenario_tree._stages[0]._tree_nodes[0] 185 186 cvar_cost_variable_name = "CVAR_COST_" + root_node._name 187 cvar_cost_variable = Var(name=cvar_cost_variable_name) 188 setattr(binding_instance, cvar_cost_variable_name, cvar_cost_variable) 189 cvar_cost_variable.construct() 190 191 cvar_eta_variable_name = "CVAR_ETA_" + root_node._name 192 cvar_eta_variable = Var(name=cvar_eta_variable_name) 193 setattr(binding_instance, cvar_eta_variable_name, cvar_eta_variable) 194 cvar_eta_variable.construct() 195 196 binding_instance.preprocess() 197 198 # ditto above for the (non-expected) cost variable. 199 for stage in scenario_tree._stages: 200 201 (cost_variable,cost_variable_index) = stage._cost_variable 202 203 print "Creating master variable and blending constraints for cost variable="+cost_variable.name+", index="+str(cost_variable_index) 204 205 for tree_node in stage._tree_nodes: 206 207 # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!) 208 209 # this is undoubtedly wasteful, in that a cost variable 210 # for each tree node is created with *all* indices. 211 new_cost_variable_name = tree_node._name + "_" + cost_variable.name 212 new_cost_variable_index = cost_variable._index 213 new_cost_variable = None 214 if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index): 215 new_cost_variable = Var(name=new_cost_variable_name) 216 else: 217 new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name) 218 new_cost_variable.construct() 219 new_cost_variable._model = binding_instance 220 setattr(binding_instance, new_cost_variable_name, new_cost_variable) 221 222 # the following is necessary, specifically to get the name - deepcopy won't reset these attributes. 223 new_cost_variable[cost_variable_index].var = new_cost_variable 224 if cost_variable_index is not None: 225 # if the variable index is None, the variable is derived from a VarValue, so the 226 # name gets updated automagically. 227 new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name 228 229 for scenario in tree_node._scenarios: 230 231 scenario_instance = scenario_instances[scenario._name] 232 scenario_cost_variable = getattr(scenario_instance, cost_variable.name) 233 new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index) 234 new_constraint = Constraint(name=new_constraint_name) 235 new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index] 236 new_constraint.add(None, (0.0, new_expr, 0.0)) 237 new_constraint._model = binding_instance 238 setattr(binding_instance, new_constraint_name, new_constraint) 239 240 # create the constraints for computing the master per-node cost variables, 241 # i.e., the current node cost and the expected cost of the child nodes. 242 # if the root, then the constraint is just the objective. 243 244 for stage in scenario_tree._stages: 245 246 (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable 247 248 for tree_node in stage._tree_nodes: 249 250 node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name 251 node_expected_cost_variable = getattr(binding_instance, node_expected_cost_variable_name) 252 253 node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name 254 node_cost_variable = getattr(binding_instance, node_cost_variable_name) 255 256 constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index] 257 258 for child_node in tree_node._children: 259 260 child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name 261 child_node_expected_cost_variable = getattr(binding_instance, child_node_expected_cost_variable_name) 262 constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable) 263 264 new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index) 265 new_constraint = Constraint(name=new_constraint_name) 266 new_constraint.add(None, (0.0, constraint_expr, 0.0)) 267 new_constraint._model = binding_instance 268 setattr(binding_instance, new_constraint_name, new_constraint) 269 270 if tree_node._parent is None: 271 272 an_instance = scenario_instances[scenario_instances.keys()[0]] 273 an_objective = an_instance.active_components(Objective) 274 opt_sense = an_objective[an_objective.keys()[0]].sense 275 276 opt_expression = node_expected_cost_variable 277 278 if generate_weighted_cvar is True: 279 cvar_cost_variable_name = "CVAR_COST_" + tree_node._name 280 cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name) 281 opt_expression += cvar_weight * cvar_cost_variable 282 283 new_objective = Objective(name="MASTER", sense=opt_sense) 284 new_objective._data[None].expr = node_expected_cost_variable 285 setattr(binding_instance, "MASTER", new_objective) 286 287 # CVaR requires the addition of a variable per scenario to represent the cost excess, 288 # and a constraint to compute the cost excess relative to eta. we also replicate (following 289 # what we do for node cost variables) an eta variable for each scenario instance, and 290 # require equality with the master eta variable via constraints. 291 if generate_weighted_cvar is True: 292 293 root_node = scenario_tree._stages[0]._tree_nodes[0] 294 295 master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name 296 master_cvar_eta_variable = getattr(binding_instance, master_cvar_eta_variable_name) 297 298 for scenario_name in scenario_instances.keys(): 299 scenario_instance = scenario_instances[scenario_name] 300 301 # unique names are required because the presolve isn't 302 # aware of the "owning" models for variables. 303 cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name 304 cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals) 305 setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable) 306 cvar_excess_variable.construct() 307 308 cvar_eta_variable_name = "CVAR_ETA" 309 cvar_eta_variable = Var(name=cvar_eta_variable_name) 310 setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable) 311 cvar_eta_variable.construct() 312 313 compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS" 314 compute_excess_constraint = Constraint(name=compute_excess_constraint_name) 315 compute_excess_expression = cvar_excess_variable 316 for node in scenario_tree._scenario_map[scenario_name]._node_list: 317 (cost_variable, cost_variable_idx) = node._stage._cost_variable 318 compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx] 319 compute_excess_expression += cvar_eta_variable 320 compute_excess_constraint.add(None, (0.0, compute_excess_expression, None)) 321 compute_excess_constraint._model = scenario_instance 322 setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint) 323 324 eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name 325 eta_equality_constraint = Constraint(name=eta_equality_constraint_name) 326 eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable 327 eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0)) 328 eta_equality_constraint._model = binding_instance 329 setattr(binding_instance, eta_equality_constraint_name, eta_equality_constraint) 330 331 # add the constraint to compute the master CVaR variable value. iterate 332 # over scenario instances to create the expected excess component first. 333 cvar_cost_variable_name = "CVAR_COST_" + root_node._name 334 cvar_cost_variable = getattr(binding_instance, cvar_cost_variable_name) 335 cvar_eta_variable_name = "CVAR_ETA_" + root_node._name 336 cvar_eta_variable = getattr(binding_instance, cvar_eta_variable_name) 337 338 cvar_cost_expression = cvar_cost_variable - cvar_eta_variable 339 340 for scenario_name in scenario_instances.keys(): 341 scenario_instance = scenario_instances[scenario_name] 342 scenario_probability = scenario_tree._scenario_map[scenario_name]._probability 343 344 scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name 345 scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name) 346 347 cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha) 348 349 compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST" 350 compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name) 351 compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0)) 352 compute_cvar_cost_constraint._model = binding_instance 353 setattr(binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint) 354 355 # after mucking with instances, presolve to collect terms required prior to output. 356 # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios. 357 for scenario_name in scenario_instances.keys(): 358 scenario_instance = scenario_instances[scenario_name] 359 scenario_instance.preprocess() 360 361 binding_instance.preprocess() 362 363 return binding_instance 364 365 # 366 # write the EF binding instance and all sub-instances. currently only outputs the CPLEX LP file format. 367 # 368 369 def write_ef(binding_instance, scenario_instances, output_filename): 370 371 # create the output file. 372 problem_writer = cpxlp.ProblemWriter_cpxlp() 373 output_file = open(output_filename,"w") 374 375 problem_writer._output_prefixes = True # we always want prefixes 376 377 ################################################################################################ 378 #### WRITE THE MASTER OBJECTIVE ################################################################ 379 ################################################################################################ 380 381 # write the objective for the master binding instance. 382 problem_writer._output_objectives = True 383 problem_writer._output_constraints = False 384 problem_writer._output_variables = False 385 386 print >>output_file, "\\ Begin objective block for master" 387 problem_writer._print_model_LP(binding_instance, output_file) 388 print >>output_file, "\\ End objective block for master" 389 print >>output_file, "" 390 391 ################################################################################################ 392 #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ######################## 393 ################################################################################################ 394 395 print >>output_file, "s.t." 396 print >>output_file, "" 397 398 problem_writer._output_objectives = False 399 problem_writer._output_constraints = True 400 problem_writer._output_variables = False 401 402 print >>output_file, "\\ Begin constraint block for master" 403 problem_writer._print_model_LP(binding_instance, output_file) 404 print >>output_file, "\\ End constraint block for master", 405 print >>output_file, "" 406 407 for scenario_name in scenario_instances.keys(): 408 instance = scenario_instances[scenario_name] 409 print >>output_file, "\\ Begin constraint block for scenario",scenario_name 410 problem_writer._print_model_LP(instance, output_file) 411 print >>output_file, "\\ End constraint block for scenario",scenario_name 412 print >>output_file, "" 413 414 ################################################################################################ 415 #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################## 416 ################################################################################################ 417 418 # write the variables for the master binding instance, and then for each scenario. 419 print >>output_file, "bounds" 420 print >>output_file, "" 421 422 problem_writer._output_objectives = False 423 problem_writer._output_constraints = False 424 problem_writer._output_variables = True 425 426 # first step: write variable bounds 427 428 problem_writer._output_continuous_variables = True 429 problem_writer._output_integer_variables = False 430 problem_writer._output_binary_variables = False 431 432 print >>output_file, "\\ Begin variable bounds block for master" 433 problem_writer._print_model_LP(binding_instance, output_file) 434 print >>output_file, "\\ End variable bounds block for master" 435 print >>output_file, "" 436 437 for scenario_name in scenario_instances.keys(): 438 instance = scenario_instances[scenario_name] 439 print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name 440 problem_writer._print_model_LP(instance, output_file) 441 print >>output_file, "\\ End variable bounds block for scenario",scenario_name 442 print >>output_file, "" 443 444 # second step: write integer indicators. 445 446 problem_writer._output_continuous_variables = False 447 problem_writer._output_integer_variables = True 448 449 if integers_present(binding_instance, scenario_instances) is True: 450 451 print >>output_file, "integer" 452 print >>output_file, "" 453 454 print >>output_file, "\\ Begin integer variable block for master" 455 problem_writer._print_model_LP(binding_instance, output_file) 456 print >>output_file, "\\ End integer variable block for master" 457 print >>output_file, "" 458 459 for scenario_name in scenario_instances.keys(): 460 instance = scenario_instances[scenario_name] 461 print >>output_file, "\\ Begin integer variable block for scenario",scenario_name 462 problem_writer._print_model_LP(instance, output_file) 463 print >>output_file, "\\ End integer variable block for scenario",scenario_name 464 print >>output_file, "" 465 466 # third step: write binary indicators. 467 468 problem_writer._output_integer_variables = False 469 problem_writer._output_binary_variables = True 470 471 if binaries_present(binding_instance, scenario_instances) is True: 472 473 print >>output_file, "binary" 474 print >>output_file, "" 475 476 print >>output_file, "\\ Begin binary variable block for master" 477 problem_writer._print_model_LP(binding_instance, output_file) 478 print >>output_file, "\\ End binary variable block for master" 479 print >>output_file, "" 480 481 for scenario_name in scenario_instances.keys(): 482 instance = scenario_instances[scenario_name] 483 print >>output_file, "\\ Begin binary variable block for scenario",scenario_name 484 problem_writer._print_model_LP(instance, output_file) 485 print >>output_file, "\\ End integer binary block for scenario",scenario_name 486 print >>output_file, "" 487 488 # wrap up. 489 print >>output_file, "end" 490 491 # clean up. 492 output_file.close() 493 494 495 # 55 496 # the main extensive-form writer routine - including read of scenarios/etc. 56 def write_ef_from_scratch(model_directory, instance_directory, output_filename, \ 497 # 498 499 def write_ef_from_scratch(model_directory, instance_directory, output_filename, 500 verbose_output, 57 501 generate_weighted_cvar, cvar_weight, risk_alpha): 58 502 … … 61 505 scenario_data_directory_name = instance_directory 62 506 63 print "Initializing extensive form writer" 64 print "" 65 66 ################################################################################################ 67 #### INITIALIZATION ############################################################################ 68 ################################################################################################ 69 70 # 71 # validate cvar options, if specified. 72 # 73 if generate_weighted_cvar is True: 74 if cvar_weight <= 0.0: 75 raise RuntimeError, "Weight of CVaR term must be >= 0.0 - value supplied="+str(cvar_weight) 76 if (risk_alpha <= 0.0) or (risk_alpha >= 1.0): 77 raise RuntimeError, "CVaR risk alpha must be between 0 and 1, exclusive - value supplied="+str(risk_alpha) 78 79 print "Writing CVaR weighted objective" 80 print "CVaR term weight="+str(cvar_weight) 81 print "CVaR alpha="+str(risk_alpha) 507 if verbose_output is True: 508 print "Initializing extensive form writer" 82 509 print "" 83 510 84 511 # 85 512 # create and populate the core model … … 115 542 except IOError: 116 543 117 print "***ERROR: Failed to load scenario reference instance data from file="+reference_ instance_filename544 print "***ERROR: Failed to load scenario reference instance data from file="+reference_scenario_filename 118 545 return 119 546 … … 135 562 # print the input tree for validation/information purposes. 136 563 # 137 scenario_tree.pprint() 564 if verbose_output is True: 565 scenario_tree.pprint() 138 566 139 567 # … … 145 573 sys.exit(1) 146 574 else: 147 print "Scenario tree is valid!" 575 if verbose_output is True: 576 print "Scenario tree is valid!" 148 577 print "" 149 578 … … 152 581 # 153 582 154 instances = {} 583 # the construction of instances takes little overhead in terms of 584 # memory potentially lost in the garbage-collection sense (mainly 585 # only that due to parsing and instance simplification/prep-processing). 586 # to speed things along, disable garbage collection if it enabled in 587 # the first place through the instance construction process. 588 # IDEA: If this becomes too much for truly large numbers of scenarios, 589 # we could manually collect every time X instances have been created. 590 591 re_enable_gc = False 592 if gc.isenabled() is True: 593 re_enable_gc = True 594 gc.disable() 595 596 scenario_instances = {} 155 597 156 598 if scenario_tree._scenario_based_data == 1: 157 print "Scenario-based instance initialization enabled" 599 if verbose_output is True: 600 print "Scenario-based instance initialization enabled" 158 601 else: 159 print "Node-based instance initialization enabled" 602 if verbose_output is True: 603 print "Node-based instance initialization enabled" 160 604 161 605 for scenario in scenario_tree._scenarios: … … 168 612 if scenario_tree._scenario_based_data == 1: 169 613 scenario_data_filename = scenario_data_directory_name + os.sep + scenario._name + ".dat" 170 # print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename171 614 scenario_instance = master_scenario_model.create(scenario_data_filename) 172 615 else: … … 176 619 while current_node is not None: 177 620 node_data_filename = scenario_data_directory_name + os.sep + current_node._name + ".dat" 178 # print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename179 621 scenario_data.add(node_data_filename) 180 622 current_node = current_node._parent … … 191 633 192 634 scenario_instance.preprocess() 193 instances[scenario._name] = scenario_instance 194 195 print "" 196 197 ################################################################################################ 198 #### CREATE THE MASTER / BINDING INSTANCE ###################################################### 199 ################################################################################################ 200 201 master_binding_instance = Model() 202 master_binding_instance.name = "MASTER" 203 204 # walk the scenario tree - create variables representing the common values for all scenarios 205 # associated with that node. the constraints will be created later. also create expected-cost 206 # variables for each node, to be computed via constraints/objectives defined in a subsequent pass. 207 # master variables are created for all nodes but those in the last stage. expected cost variables 208 # are, for no particularly good reason other than easy coding, created for nodes in all stages. 209 print "Creating variables for master binding instance" 210 211 for stage in scenario_tree._stages: 212 213 for (stage_variable, index_template, stage_variable_indices) in stage._variables: 214 215 print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=",stage_variable_indices 216 217 for tree_node in stage._tree_nodes: 218 219 if stage != scenario_tree._stages[-1]: 220 221 master_variable_name = tree_node._name + "_" + stage_variable.name 222 223 # because there may be a single stage variable and multiple indices, check 224 # for the existence of the variable at this node - if you don't, you'll 225 # inadvertently over-write what was there previously! 226 master_variable = None 227 try: 228 master_variable = getattr(master_binding_instance, master_variable_name) 229 except: 230 new_master_variable_index = stage_variable._index 231 new_master_variable = None 232 if (len(new_master_variable_index) is 1) and (None in new_master_variable_index): 233 new_master_variable = Var(name=stage_variable.name) 234 else: 235 new_master_variable = Var(new_master_variable_index, name=stage_variable.name) 236 new_master_variable.construct() 237 new_master_variable._model = master_binding_instance 238 setattr(master_binding_instance, master_variable_name, new_master_variable) 239 240 # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS! 241 242 master_variable = new_master_variable 243 244 for index in stage_variable_indices: 245 246 is_used = True # until proven otherwise 247 for scenario in tree_node._scenarios: 248 instance = instances[scenario._name] 249 if getattr(instance,stage_variable.name)[index].status == VarStatus.unused: 250 is_used = False 251 252 is_fixed = False # until proven otherwise 253 for scenario in tree_node._scenarios: 254 instance = instances[scenario._name] 255 if getattr(instance,stage_variable.name)[index].fixed is True: 256 is_fixed = True 257 258 if (is_used is True) and (is_fixed is False): 259 260 # the following is necessary, specifically to get the name - deepcopy won't reset these attributes. 261 # and because presolve/simplification is name-based, the names *have* to be different. 262 master_variable[index].var = master_variable 263 master_variable[index].name = tree_node._name + "_" + master_variable[index].name 264 265 for scenario in tree_node._scenarios: 266 267 scenario_instance = instances[scenario._name] 268 scenario_variable = getattr(scenario_instance, stage_variable.name) 269 new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index) 270 new_constraint = Constraint(name=new_constraint_name) 271 new_expr = master_variable[index] - scenario_variable[index] 272 new_constraint.add(None, (0.0, new_expr, 0.0)) 273 new_constraint._model = master_binding_instance 274 setattr(master_binding_instance, new_constraint_name, new_constraint) 275 276 # create a variable to represent the expected cost at this node - 277 # the constraint to compute this comes later. 278 expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name 279 expected_cost_variable = Var(name=expected_cost_variable_name) 280 expected_cost_variable._model = master_binding_instance 281 setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable) 282 283 # if we're generating the weighted CVaR objective term, create the corresponding variable and 284 # the master CVaR eta variable. 285 if generate_weighted_cvar is True: 286 root_node = scenario_tree._stages[0]._tree_nodes[0] 287 288 cvar_cost_variable_name = "CVAR_COST_" + root_node._name 289 cvar_cost_variable = Var(name=cvar_cost_variable_name) 290 setattr(master_binding_instance, cvar_cost_variable_name, cvar_cost_variable) 291 cvar_cost_variable.construct() 292 293 cvar_eta_variable_name = "CVAR_ETA_" + root_node._name 294 cvar_eta_variable = Var(name=cvar_eta_variable_name) 295 setattr(master_binding_instance, cvar_eta_variable_name, cvar_eta_variable) 296 cvar_eta_variable.construct() 297 298 master_binding_instance.preprocess() 299 300 # ditto above for the (non-expected) cost variable. 301 for stage in scenario_tree._stages: 302 303 (cost_variable,cost_variable_index) = stage._cost_variable 304 305 print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index 306 307 for tree_node in stage._tree_nodes: 308 309 # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!) 310 311 # this is undoubtedly wasteful, in that a cost variable 312 # for each tree node is created with *all* indices. 313 new_cost_variable_name = tree_node._name + "_" + cost_variable.name 314 new_cost_variable_index = cost_variable._index 315 new_cost_variable = None 316 if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index): 317 new_cost_variable = Var(name=new_cost_variable_name) 318 else: 319 new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name) 320 new_cost_variable.construct() 321 new_cost_variable._model = master_binding_instance 322 setattr(master_binding_instance, new_cost_variable_name, new_cost_variable) 323 324 # the following is necessary, specifically to get the name - deepcopy won't reset these attributes. 325 new_cost_variable[cost_variable_index].var = new_cost_variable 326 if cost_variable_index is not None: 327 # if the variable index is None, the variable is derived from a VarValue, so the 328 # name gets updated automagically. 329 new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name 330 331 for scenario in tree_node._scenarios: 332 333 scenario_instance = instances[scenario._name] 334 scenario_cost_variable = getattr(scenario_instance, cost_variable.name) 335 new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index) 336 new_constraint = Constraint(name=new_constraint_name) 337 new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index] 338 new_constraint.add(None, (0.0, new_expr, 0.0)) 339 new_constraint._model = master_binding_instance 340 setattr(master_binding_instance, new_constraint_name, new_constraint) 341 342 # create the constraints for computing the master per-node cost variables, 343 # i.e., the current node cost and the expected cost of the child nodes. 344 # if the root, then the constraint is just the objective. 345 346 for stage in scenario_tree._stages: 347 348 (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable 349 350 for tree_node in stage._tree_nodes: 351 352 node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name 353 node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name) 354 355 node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name 356 node_cost_variable = getattr(master_binding_instance, node_cost_variable_name) 357 358 constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index] 359 360 for child_node in tree_node._children: 361 362 child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name 363 child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name) 364 constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable) 365 366 new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index) 367 new_constraint = Constraint(name=new_constraint_name) 368 new_constraint.add(None, (0.0, constraint_expr, 0.0)) 369 new_constraint._model = master_binding_instance 370 setattr(master_binding_instance, new_constraint_name, new_constraint) 371 372 if tree_node._parent is None: 373 374 an_instance = instances[instances.keys()[0]] 375 an_objective = an_instance.active_components(Objective) 376 opt_sense = an_objective[an_objective.keys()[0]].sense 377 378 opt_expression = node_expected_cost_variable 379 380 if generate_weighted_cvar is True: 381 cvar_cost_variable_name = "CVAR_COST_" + tree_node._name 382 cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name) 383 opt_expression += cvar_weight * cvar_cost_variable 384 385 new_objective = Objective(name="MASTER", sense=opt_sense) 386 new_objective._data[None].expr = opt_expression 387 setattr(master_binding_instance, "MASTER", new_objective) 388 389 # CVaR requires the addition of a variable per scenario to represent the cost excess, 390 # and a constraint to compute the cost excess relative to eta. we also replicate (following 391 # what we do for node cost variables) an eta variable for each scenario instance, and 392 # require equality with the master eta variable via constraints. 393 if generate_weighted_cvar is True: 394 395 root_node = scenario_tree._stages[0]._tree_nodes[0] 396 397 master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name 398 master_cvar_eta_variable = getattr(master_binding_instance, master_cvar_eta_variable_name) 399 400 for scenario_name in instances.keys(): 401 scenario_instance = instances[scenario_name] 402 403 # unique names are required because the presolve isn't 404 # aware of the "owning" models for variables. 405 cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name 406 cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals) 407 setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable) 408 cvar_excess_variable.construct() 409 410 cvar_eta_variable_name = "CVAR_ETA" 411 cvar_eta_variable = Var(name=cvar_eta_variable_name) 412 setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable) 413 cvar_eta_variable.construct() 414 415 compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS" 416 compute_excess_constraint = Constraint(name=compute_excess_constraint_name) 417 compute_excess_expression = cvar_excess_variable 418 for node in scenario_tree._scenario_map[scenario_name]._node_list: 419 (cost_variable, cost_variable_idx) = node._stage._cost_variable 420 compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx] 421 compute_excess_expression += cvar_eta_variable 422 compute_excess_constraint.add(None, (0.0, compute_excess_expression, None)) 423 compute_excess_constraint._model = scenario_instance 424 setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint) 425 426 eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name 427 eta_equality_constraint = Constraint(name=eta_equality_constraint_name) 428 eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable 429 eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0)) 430 eta_equality_constraint._model = master_binding_instance 431 setattr(master_binding_instance, eta_equality_constraint_name, eta_equality_constraint) 432 433 # add the constraint to compute the master CVaR variable value. iterate 434 # over scenario instances to create the expected excess component first. 435 cvar_cost_variable_name = "CVAR_COST_" + root_node._name 436 cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name) 437 cvar_eta_variable_name = "CVAR_ETA_" + root_node._name 438 cvar_eta_variable = getattr(master_binding_instance, cvar_eta_variable_name) 439 440 cvar_cost_expression = cvar_cost_variable - cvar_eta_variable 441 442 for scenario_name in instances.keys(): 443 scenario_instance = instances[scenario_name] 444 scenario_probability = scenario_tree._scenario_map[scenario_name]._probability 445 446 scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name 447 scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name) 448 449 cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha) 450 451 compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST" 452 compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name) 453 compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0)) 454 compute_cvar_cost_constraint._model = master_binding_instance 455 setattr(master_binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint) 456 457 # after mucking with instances, presolve to collect terms required prior to output. 458 # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios. 459 for scenario_name in instances.keys(): 460 scenario_instance = instances[scenario_name] 461 scenario_instance.preprocess() 462 463 master_binding_instance.preprocess() 464 465 ################################################################################################ 466 #### WRITE THE COMPOSITE MODEL ################################################################# 467 ################################################################################################ 635 scenario_instances[scenario._name] = scenario_instance 636 637 if re_enable_gc is True: 638 gc.enable() 639 640 print "" 641 642 binding_instance = create_ef_instance(scenario_tree, scenario_instances, 643 verbose_output = verbose_output, 644 generate_weighted_cvar = generate_weighted_cvar, 645 cvar_weight = cvar_weight, 646 risk_alpha = risk_alpha) 468 647 469 648 print "" 470 649 print "Starting to write extensive form" 471 650 472 # create the output file. 473 problem_writer = cpxlp.ProblemWriter_cpxlp() 474 output_file = open(output_filename,"w") 475 476 problem_writer._output_prefixes = True # we always want prefixes 477 478 ################################################################################################ 479 #### WRITE THE MASTER OBJECTIVE ################################################################ 480 ################################################################################################ 481 482 # write the objective for the master binding instance. 483 problem_writer._output_objectives = True 484 problem_writer._output_constraints = False 485 problem_writer._output_variables = False 486 487 print >>output_file, "\\ Begin objective block for master" 488 problem_writer._print_model_LP(master_binding_instance, output_file) 489 print >>output_file, "\\ End objective block for master" 490 print >>output_file, "" 491 492 ################################################################################################ 493 #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ######################## 494 ################################################################################################ 495 496 print >>output_file, "s.t." 497 print >>output_file, "" 498 499 problem_writer._output_objectives = False 500 problem_writer._output_constraints = True 501 problem_writer._output_variables = False 502 503 print >>output_file, "\\ Begin constraint block for master" 504 problem_writer._print_model_LP(master_binding_instance, output_file) 505 print >>output_file, "\\ End constraint block for master", 506 print >>output_file, "" 507 508 for scenario_name in instances.keys(): 509 instance = instances[scenario_name] 510 print >>output_file, "\\ Begin constraint block for scenario",scenario_name 511 problem_writer._print_model_LP(instance, output_file) 512 print >>output_file, "\\ End constraint block for scenario",scenario_name 513 print >>output_file, "" 514 515 ################################################################################################ 516 #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################## 517 ################################################################################################ 518 519 # write the variables for the master binding instance, and then for each scenario. 520 print >>output_file, "bounds" 521 print >>output_file, "" 522 523 problem_writer._output_objectives = False 524 problem_writer._output_constraints = False 525 problem_writer._output_variables = True 526 527 # first step: write variable bounds 528 529 problem_writer._output_continuous_variables = True 530 problem_writer._output_integer_variables = False 531 problem_writer._output_binary_variables = False 532 533 print >>output_file, "\\ Begin variable bounds block for master" 534 problem_writer._print_model_LP(master_binding_instance, output_file) 535 print >>output_file, "\\ End variable bounds block for master" 536 print >>output_file, "" 537 538 for scenario_name in instances.keys(): 539 instance = instances[scenario_name] 540 print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name 541 problem_writer._print_model_LP(instance, output_file) 542 print >>output_file, "\\ End variable bounds block for scenario",scenario_name 543 print >>output_file, "" 544 545 # second step: write integer indicators. 546 547 problem_writer._output_continuous_variables = False 548 problem_writer._output_integer_variables = True 549 550 if integers_present(master_binding_instance, instances) is True: 551 552 print >>output_file, "integer" 553 print >>output_file, "" 554 555 print >>output_file, "\\ Begin integer variable block for master" 556 problem_writer._print_model_LP(master_binding_instance, output_file) 557 print >>output_file, "\\ End integer variable block for master" 558 print >>output_file, "" 559 560 for scenario_name in instances.keys(): 561 instance = instances[scenario_name] 562 print >>output_file, "\\ Begin integer variable block for scenario",scenario_name 563 problem_writer._print_model_LP(instance, output_file) 564 print >>output_file, "\\ End integer variable block for scenario",scenario_name 565 print >>output_file, "" 566 567 # third step: write binary indicators. 568 569 problem_writer._output_integer_variables = False 570 problem_writer._output_binary_variables = True 571 572 if binaries_present(master_binding_instance, instances) is True: 573 574 print >>output_file, "binary" 575 print >>output_file, "" 576 577 print >>output_file, "\\ Begin binary variable block for master" 578 problem_writer._print_model_LP(master_binding_instance, output_file) 579 print >>output_file, "\\ End binary variable block for master" 580 print >>output_file, "" 581 582 for scenario_name in instances.keys(): 583 instance = instances[scenario_name] 584 print >>output_file, "\\ Begin binary variable block for scenario",scenario_name 585 problem_writer._print_model_LP(instance, output_file) 586 print >>output_file, "\\ End integer binary block for scenario",scenario_name 587 print >>output_file, "" 588 589 # wrap up. 590 print >>output_file, "end" 591 592 # clean up. 593 output_file.close() 651 write_ef(binding_instance, scenario_instances, output_filename) 594 652 595 653 print "" … … 605 663 print "" 606 664 607 def write_ef(scenario_tree,instances, output_filename):665 def create_and_write_ef(scenario_tree, scenario_instances, output_filename): 608 666 609 667 start_time = time.time() 610 668 611 ################################################################################################ 612 #### CREATE THE MASTER / BINDING INSTANCE ###################################################### 613 ################################################################################################ 614 615 master_binding_instance = Model() 616 master_binding_instance.name = "MASTER" 617 618 # walk the scenario tree - create variables representing the common values for all scenarios 619 # associated with that node. the constraints will be created later. also create expected-cost 620 # variables for each node, to be computed via constraints/objectives defined in a subsequent pass. 621 # master variables are created for all nodes but those in the last stage. expected cost variables 622 # are, for no particularly good reason other than easy coding, created for nodes in all stages. 623 print "Creating variables for master binding instance" 624 625 for stage in scenario_tree._stages: 626 627 for (stage_variable, index_template, stage_variable_indices) in stage._variables: 628 629 print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=", index_template 630 631 for tree_node in stage._tree_nodes: 632 633 if stage != scenario_tree._stages[-1]: 634 635 master_variable_name = tree_node._name + "_" + stage_variable.name 636 637 # because there may be a single stage variable and multiple indices, check 638 # for the existence of the variable at this node - if you don't, you'll 639 # inadvertently over-write what was there previously! 640 master_variable = None 641 try: 642 master_variable = getattr(master_binding_instance, master_variable_name) 643 except: 644 new_master_variable_index = stage_variable._index 645 new_master_variable = None 646 if (len(new_master_variable_index) is 1) and (None in new_master_variable_index): 647 new_master_variable = Var(name=stage_variable.name) 648 else: 649 new_master_variable = Var(new_master_variable_index, name=stage_variable.name) 650 new_master_variable.construct() 651 new_master_variable._model = master_binding_instance 652 setattr(master_binding_instance, master_variable_name, new_master_variable) 653 654 # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS! 655 656 master_variable = new_master_variable 657 658 for index in stage_variable_indices: 659 660 is_used = True # until proven otherwise 661 for scenario in tree_node._scenarios: 662 instance = instances[scenario._name] 663 if getattr(instance,stage_variable.name)[index].status == VarStatus.unused: 664 is_used = False 665 666 is_fixed = False # until proven otherwise 667 for scenario in tree_node._scenarios: 668 instance = instances[scenario._name] 669 if getattr(instance,stage_variable.name)[index].fixed is True: 670 is_fixed = True 671 672 if (is_used is True) and (is_fixed is False): 673 674 # the following is necessary, specifically to get the name - deepcopy won't reset these attributes. 675 # and because presolve/simplification is name-based, the names *have* to be different. 676 master_variable[index].var = master_variable 677 master_variable[index].name = tree_node._name + "_" + master_variable[index].name 678 679 for scenario in tree_node._scenarios: 680 681 scenario_instance = instances[scenario._name] 682 scenario_variable = getattr(scenario_instance, stage_variable.name) 683 new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index) 684 new_constraint = Constraint(name=new_constraint_name) 685 new_expr = master_variable[index] - scenario_variable[index] 686 new_constraint.add(None, (0.0, new_expr, 0.0)) 687 new_constraint._model = master_binding_instance 688 setattr(master_binding_instance, new_constraint_name, new_constraint) 689 690 # create a variable to represent the expected cost at this node - 691 # the constraint to compute this comes later. 692 expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name 693 expected_cost_variable = Var(name=expected_cost_variable_name) 694 expected_cost_variable._model = master_binding_instance 695 setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable) 696 697 master_binding_instance.preprocess() 698 699 # ditto above for the (non-expected) cost variable. 700 for stage in scenario_tree._stages: 701 702 (cost_variable,cost_variable_index) = stage._cost_variable 703 704 print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index 705 706 for tree_node in stage._tree_nodes: 707 708 new_cost_variable_name = tree_node._name + "_" + cost_variable.name 709 710 # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!) 711 712 # this is undoubtedly wasteful, in that a cost variable 713 # for each tree node is created with *all* indices. 714 new_cost_variable_name = tree_node._name + "_" + cost_variable.name 715 new_cost_variable_index = cost_variable._index 716 new_cost_variable = None 717 if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index): 718 new_cost_variable = Var(name=new_cost_variable_name) 719 else: 720 new_cost_variable = Var(new_cost_variable_index, new_cost_variable_name) 721 new_cost_variable.construct() 722 new_cost_variable._model = master_binding_instance 723 setattr(master_binding_instance, new_cost_variable_name, new_cost_variable) 724 725 # the following is necessary, specifically to get the name - deepcopy won't reset these attributes. 726 new_cost_variable[cost_variable_index].var = new_cost_variable 727 if cost_variable_index is not None: 728 # if the variable index is None, the variable is derived from a VarValue, so the 729 # name gets updated automagically. 730 new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name 731 732 for scenario in tree_node._scenarios: 733 734 scenario_instance = instances[scenario._name] 735 scenario_cost_variable = getattr(scenario_instance, cost_variable.name) 736 new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index) 737 new_constraint = Constraint(name=new_constraint_name) 738 new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index] 739 new_constraint.add(None, (0.0, new_expr, 0.0)) 740 new_constraint._model = master_binding_instance 741 setattr(master_binding_instance, new_constraint_name, new_constraint) 742 743 # create the constraints for computing the master per-node cost variables, 744 # i.e., the current node cost and the expected cost of the child nodes. 745 # if the root, then the constraint is just the objective. 746 747 for stage in scenario_tree._stages: 748 749 (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable 750 751 for tree_node in stage._tree_nodes: 752 753 node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name 754 node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name) 755 756 node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name 757 node_cost_variable = getattr(master_binding_instance, node_cost_variable_name) 758 759 constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index] 760 761 for child_node in tree_node._children: 762 763 child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name 764 child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name) 765 constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable) 766 767 new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index) 768 new_constraint = Constraint(name=new_constraint_name) 769 new_constraint.add(None, (0.0, constraint_expr, 0.0)) 770 new_constraint._model = master_binding_instance 771 setattr(master_binding_instance, new_constraint_name, new_constraint) 772 773 if tree_node._parent is None: 774 775 an_instance = instances[instances.keys()[0]] 776 an_objective = an_instance.active_components(Objective) 777 opt_sense = an_objective[an_objective.keys()[0]].sense 778 779 new_objective = Objective(name="MASTER", sense=opt_sense) 780 new_objective._data[None].expr = node_expected_cost_variable 781 setattr(master_binding_instance, "MASTER", new_objective) 782 783 master_binding_instance.preprocess() 784 785 ################################################################################################ 786 #### WRITE THE COMPOSITE MODEL ################################################################# 787 ################################################################################################ 669 binding_instance = create_ef_instance(scenario_tree, scenario_instances) 788 670 789 671 print "" 790 672 print "Starting to write extensive form" 791 673 792 # create the output file. 793 problem_writer = cpxlp.ProblemWriter_cpxlp() 794 output_file = open(output_filename,"w") 795 796 problem_writer._output_prefixes = True # we always want prefixes 797 798 ################################################################################################ 799 #### WRITE THE MASTER OBJECTIVE ################################################################ 800 ################################################################################################ 801 802 # write the objective for the master binding instance. 803 problem_writer._output_objectives = True 804 problem_writer._output_constraints = False 805 problem_writer._output_variables = False 806 807 print >>output_file, "\\ Begin objective block for master" 808 problem_writer._print_model_LP(master_binding_instance, output_file) 809 print >>output_file, "\\ End objective block for master" 810 print >>output_file, "" 811 812 ################################################################################################ 813 #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ######################## 814 ################################################################################################ 815 816 print >>output_file, "s.t." 817 print >>output_file, "" 818 819 problem_writer._output_objectives = False 820 problem_writer._output_constraints = True 821 problem_writer._output_variables = False 822 823 print >>output_file, "\\ Begin constraint block for master" 824 problem_writer._print_model_LP(master_binding_instance, output_file) 825 print >>output_file, "\\ End constraint block for master", 826 print >>output_file, "" 827 828 for scenario_name in instances.keys(): 829 instance = instances[scenario_name] 830 print >>output_file, "\\ Begin constraint block for scenario",scenario_name 831 problem_writer._print_model_LP(instance, output_file) 832 print >>output_file, "\\ End constraint block for scenario",scenario_name 833 print >>output_file, "" 834 835 ################################################################################################ 836 #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################## 837 ################################################################################################ 838 839 # write the variables for the master binding instance, and then for each scenario. 840 print >>output_file, "bounds" 841 print >>output_file, "" 842 843 problem_writer._output_objectives = False 844 problem_writer._output_constraints = False 845 problem_writer._output_variables = True 846 847 # first step: write variable bounds 848 849 problem_writer._output_continuous_variables = True 850 problem_writer._output_integer_variables = False 851 problem_writer._output_binary_variables = False 852 853 print >>output_file, "\\ Begin variable bounds block for master" 854 problem_writer._print_model_LP(master_binding_instance, output_file) 855 print >>output_file, "\\ End variable bounds block for master" 856 print >>output_file, "" 857 858 for scenario_name in instances.keys(): 859 instance = instances[scenario_name] 860 print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name 861 problem_writer._print_model_LP(instance, output_file) 862 print >>output_file, "\\ End variable bounds block for scenario",scenario_name 863 print >>output_file, "" 864 865 # second step: write integer indicators. 866 867 problem_writer._output_continuous_variables = False 868 problem_writer._output_integer_variables = True 869 870 if integers_present(master_binding_instance, instances) is True: 871 872 print >>output_file, "integer" 873 print >>output_file, "" 874 875 print >>output_file, "\\ Begin integer variable block for master" 876 problem_writer._print_model_LP(master_binding_instance, output_file) 877 print >>output_file, "\\ End integer variable block for master" 878 print >>output_file, "" 879 880 for scenario_name in instances.keys(): 881 instance = instances[scenario_name] 882 print >>output_file, "\\ Begin integer variable block for scenario",scenario_name 883 problem_writer._print_model_LP(instance, output_file) 884 print >>output_file, "\\ End integer variable block for scenario",scenario_name 885 print >>output_file, "" 886 887 # third step: write binary indicators. 888 889 problem_writer._output_integer_variables = False 890 problem_writer._output_binary_variables = True 891 892 if binaries_present(master_binding_instance, instances) is True: 893 894 print >>output_file, "binary" 895 print >>output_file, "" 896 897 print >>output_file, "\\ Begin binary variable block for master" 898 problem_writer._print_model_LP(master_binding_instance, output_file) 899 print >>output_file, "\\ End binary variable block for master" 900 print >>output_file, "" 901 902 for scenario_name in instances.keys(): 903 instance = instances[scenario_name] 904 print >>output_file, "\\ Begin binary variable block for scenario",scenario_name 905 problem_writer._print_model_LP(instance, output_file) 906 print >>output_file, "\\ End integer binary block for scenario",scenario_name 907 print >>output_file, "" 908 909 # wrap up. 910 print >>output_file, "end" 911 912 # clean up. 913 output_file.close() 674 write_ef(binding_instance, scenario_instances, output_filename) 914 675 915 676 print "" -
coopr.pysp/stable/2.3/coopr/pysp/ef_writer_script.py
r2317 r2434 24 24 25 25 # 26 # Setup command-line options26 # utility method to construct an option parser for ef writer arguments 27 27 # 28 28 29 parser = OptionParser() 30 parser.add_option("--model-directory", 31 help="The directory in which all model (reference and scenario) definitions are stored. Default is \".\".", 32 action="store", 33 dest="model_directory", 34 type="string", 35 default=".") 36 parser.add_option("--instance-directory", 37 help="The directory in which all instance (reference and scenario) definitions are stored. Default is \".\".", 38 action="store", 39 dest="instance_directory", 40 type="string", 41 default=".") 42 parser.add_option("--generate-weighted-cvar", 43 help="Add a weighted CVaR term to the primary objective", 44 action="store_true", 45 dest="generate_weighted_cvar", 46 default=False) 47 parser.add_option("--cvar-weight", 48 help="The weight associated with the CVaR term in the risk-weighted objective formulation.", 49 action="store", 50 dest="cvar_weight", 51 type="float", 52 default=1.0) 53 parser.add_option("--risk-alpha", 54 help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics.", 55 action="store", 56 dest="risk_alpha", 57 type="float", 58 default=0.95) 59 parser.add_option("--output-file", 60 help="Specify the name of the extensive form output file", 61 action="store", 62 dest="output_file", 63 type="string", 64 default="efout.lp") 65 parser.add_option("--profile", 66 help="Enable profiling of Python code. The value of this option is the number of functions that are summarized.", 67 action="store", 68 dest="profile", 69 default=0) 70 parser.usage="efwriter [options]" 29 def construct_ef_writer_options_parser(usage_string): 71 30 31 parser = OptionParser() 32 parser.add_option("--verbose", 33 help="Generate verbose output, beyond the usual status output. Default is False.", 34 action="store_true", 35 dest="verbose", 36 default=False) 37 parser.add_option("--model-directory", 38 help="The directory in which all model (reference and scenario) definitions are stored. Default is \".\".", 39 action="store", 40 dest="model_directory", 41 type="string", 42 default=".") 43 parser.add_option("--instance-directory", 44 help="The directory in which all instance (reference and scenario) definitions are stored. Default is \".\".", 45 action="store", 46 dest="instance_directory", 47 type="string", 48 default=".") 49 parser.add_option("--generate-weighted-cvar", 50 help="Add a weighted CVaR term to the primary objective", 51 action="store_true", 52 dest="generate_weighted_cvar", 53 default=False) 54 parser.add_option("--cvar-weight", 55 help="The weight associated with the CVaR term in the risk-weighted objective formulation.", 56 action="store", 57 dest="cvar_weight", 58 type="float", 59 default=1.0) 60 parser.add_option("--risk-alpha", 61 help="The probability threshold associated with cvar (or any future) risk-oriented performance metrics.", 62 action="store", 63 dest="risk_alpha", 64 type="float", 65 default=0.95) 66 parser.add_option("--output-file", 67 help="Specify the name of the extensive form output file", 68 action="store", 69 dest="output_file", 70 type="string", 71 default="efout.lp") 72 parser.add_option("--profile", 73 help="Enable profiling of Python code. The value of this option is the number of functions that are summarized.", 74 action="store", 75 dest="profile", 76 default=0) 77 parser.add_option("--disable-gc", 78 help="Disable the python garbage collecter. Default is False.", 79 action="store_true", 80 dest="disable_gc", 81 default=False) 82 parser.usage=usage_string 83 84 return parser 85 72 86 def run_ef_writer(options, args): 73 87 … … 84 98 risk_alpha = options.risk_alpha 85 99 86 write_ef_from_scratch(os.path.expanduser(options.model_directory), os.path.expanduser(options.instance_directory), os.path.expanduser(options.output_file), \ 100 write_ef_from_scratch(os.path.expanduser(options.model_directory), 101 os.path.expanduser(options.instance_directory), 102 os.path.expanduser(options.output_file), 103 options.verbose, 87 104 generate_weighted_cvar, cvar_weight, risk_alpha) 88 89 return90 105 91 106 def run(args=None): … … 96 111 # 97 112 98 # for a one-pass execution, garbage collection doesn't make99 # much sense - so I'm disabling it. Because: It drops the run-time100 # by a factor of 3-4 on bigger instances.101 gc.disable()102 103 113 # 104 114 # Parse command-line options. 105 115 # 106 116 try: 107 (options, args) = parser.parse_args(args=args) 117 options_parser = construct_ef_writer_options_parser("efwriter [options]") 118 (options, args) = options_parser.parse_args(args=args) 108 119 except SystemExit: 109 120 # the parser throws a system exit if "-h" is specified - catch 110 121 # it to exit gracefully. 111 122 return 123 124 if options.disable_gc is True: 125 gc.disable() 126 else: 127 gc.enable() 112 128 113 129 if options.profile > 0: … … 140 156 141 157 gc.enable() 158 142 159 return ans 143 160 -
coopr.pysp/stable/2.3/coopr/pysp/ph.py
r2391 r2434 27 27 from scenariotree import * 28 28 from phutils import * 29 from phobjective import * 29 30 30 31 from pyutilib.component.core import ExtensionPoint … … 33 34 34 35 class ProgressiveHedging(object): 35 36 #37 # routine to compute linearization breakpoints uniformly between the bounds and the mean.38 #39 40 # IMPT: In general, the breakpoint computation codes can return a 2-list even if the lb equals41 # the ub. This case happens quite often in real models (although typically lb=xvag=ub).42 # See the code for constructing the pieces on how this case is handled in the linearization.43 44 def compute_uniform_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):45 46 breakpoints = []47 48 # add the lower bound - the first breakpoint.49 breakpoints.append(lb)50 51 # determine the breakpoints to the left of the mean.52 left_step = (xavg - lb) / num_breakpoints_per_side53 current_x = lb54 for i in range(1,num_breakpoints_per_side+1):55 this_lb = current_x56 this_ub = current_x+left_step57 if (fabs(this_lb-lb) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):58 breakpoints.append(this_lb)59 current_x += left_step60 61 # add the mean - it's always a breakpoint. unless!62 # the lb or ub and the avg are the same.63 if (fabs(lb-xavg) > self._integer_tolerance) and (fabs(ub-xavg) > self._integer_tolerance):64 breakpoints.append(xavg)65 66 # determine the breakpoints to the right of the mean.67 right_step = (ub - xavg) / num_breakpoints_per_side68 current_x = xavg69 for i in range(1,num_breakpoints_per_side+1):70 this_lb = current_x71 this_ub = current_x+right_step72 if (fabs(this_ub-xavg) > self._integer_tolerance) and (fabs(this_ub-ub) > self._integer_tolerance):73 breakpoints.append(this_ub)74 current_x += right_step75 76 # add the upper bound - the last breakpoint.77 # the upper bound should always be different than the lower bound (I say with some78 # hesitation - it really depends on what plugins are doing to modify the bounds dynamically).79 breakpoints.append(ub)80 81 return breakpoints82 83 #84 # routine to compute linearization breakpoints uniformly between the current node min/max bounds.85 #86 87 def compute_uniform_between_nodestat_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):88 89 breakpoints = []90 91 # add the lower bound - the first breakpoint.92 breakpoints.append(lb)93 94 # add the node-min - the second breakpoint. but only if it is different than the lower bound and the mean.95 if (fabs(node_min-lb) > self._integer_tolerance) and (fabs(node_min-xavg) > self._integer_tolerance):96 breakpoints.append(node_min)97 98 step = (node_max - node_min) / num_breakpoints99 current_x = node_min100 for i in range(1,num_breakpoints+1):101 this_lb = current_x102 this_ub = current_x+step103 if (fabs(this_lb-node_min) > self._integer_tolerance) and (fabs(this_lb-node_max) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):104 breakpoints.append(this_lb)105 current_x += step106 107 # add the node-max - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.108 if (fabs(node_max-ub) > self._integer_tolerance) and (fabs(node_max-xavg) > self._integer_tolerance):109 breakpoints.append(node_max)110 111 # add the upper bound - the last breakpoint.112 breakpoints.append(ub)113 114 # add the mean - it's always a breakpoint. unless! -115 # it happens to be equal to (within tolerance) the lower or upper bounds.116 # sort to insert it in the correct place.117 if (fabs(xavg - lb) > self._integer_tolerance) and (fabs(xavg - ub) > self._integer_tolerance):118 breakpoints.append(xavg)119 breakpoints.sort()120 121 return breakpoints122 123 #124 # routine to compute linearization breakpoints using "Woodruff" relaxation of the compute_uniform_between_nodestat_breakpoints.125 #126 127 def compute_uniform_between_woodruff_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):128 129 breakpoints = []130 131 # add the lower bound - the first breakpoint.132 breakpoints.append(lb)133 134 # be either three "differences" from the mean, or else "halfway to the bound", whichever is closer to the mean.135 left = max(xavg - 3.0 * (xavg - node_min), xavg - 0.5 * (xavg - lb))136 right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))137 138 # add the left bound - the second breakpoint. but only if it is different than the lower bound and the mean.139 if (fabs(left-lb) > self._integer_tolerance) and (fabs(left-xavg) > self._integer_tolerance):140 breakpoints.append(left)141 142 step = (right - left) / num_breakpoints143 current_x = left144 for i in range(1,num_breakpoints+1):145 this_lb = current_x146 this_ub = current_x+step147 if (fabs(this_lb-left) > self._integer_tolerance) and (fabs(this_lb-right) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):148 breakpoints.append(this_lb)149 current_x += step150 151 # add the right bound - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.152 if (fabs(right-ub) > self._integer_tolerance) and (fabs(right-xavg) > self._integer_tolerance):153 breakpoints.append(right)154 155 # add the upper bound - the last breakpoint.156 breakpoints.append(ub)157 158 # add the mean - it's always a breakpoint.159 # sort to insert it in the correct place.160 breakpoints.append(xavg)161 breakpoints.sort()162 163 return breakpoints164 165 #166 # routine to compute linearization breakpoints based on an exponential distribution from the mean in each direction.167 #168 169 def compute_exponential_from_mean_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):170 171 breakpoints = []172 173 # add the lower bound - the first breakpoint.174 breakpoints.append(lb)175 176 # determine the breakpoints to the left of the mean.177 left_delta = xavg - lb178 base = exp(log(left_delta) / num_breakpoints_per_side)179 current_offset = base180 for i in range(1,num_breakpoints_per_side+1):181 current_x = xavg - current_offset182 if (fabs(current_x-lb) > self._integer_tolerance) and (fabs(current_x-xavg) > self._integer_tolerance):183 breakpoints.append(current_x)184 current_offset *= base185 186 # add the mean - it's always a breakpoint.187 breakpoints.append(xavg)188 189 # determine the breakpoints to the right of the mean.190 right_delta = ub - xavg191 base = exp(log(right_delta) / num_breakpoints_per_side)192 current_offset = base193 for i in range(1,num_breakpoints_per_side+1):194 current_x = xavg + current_offset195 if (fabs(current_x-xavg) > self._integer_tolerance) and (fabs(current_x-ub) > self._integer_tolerance):196 breakpoints.append(current_x)197 current_offset *= base198 199 # add the upper bound - the last breakpoint.200 breakpoints.append(ub)201 202 return breakpoints203 36 204 37 # … … 379 212 return (num_fixed_discrete_vars, num_fixed_continuous_vars) 380 213 381 # a utility to create piece-wise linear constraint expressions for a given variable, for382 # use in constructing the augmented (penalized) PH objective. lb and ub are the bounds383 # on this piece, variable is the actual instance variable, and average is the instance384 # parameter specifying the average of this variable across instances sharing passing385 # through a common tree node. lb and ub are floats.386 # IMPT: There are cases where lb=ub, in which case the slope is 0 and the intercept387 # is simply the penalty at the lower(or upper) bound.388 def _create_piecewise_constraint_expression(self, lb, ub, instance_variable, variable_average, quad_variable):389 390 penalty_at_lb = (lb - variable_average()) * (lb - variable_average())391 penalty_at_ub = (ub - variable_average()) * (ub - variable_average())392 slope = None393 if fabs(ub-lb) > self._integer_tolerance:394 slope = (penalty_at_ub - penalty_at_lb) / (ub - lb)395 else:396 slope = 0.0397 intercept = penalty_at_lb - slope * lb398 expression = (0.0, quad_variable - slope * instance_variable - intercept, None)399 400 return expression401 402 214 # when the quadratic penalty terms are approximated via piecewise linear segments, 403 215 # we end up (necessarily) "littering" the scenario instances with extra constraints. … … 424 236 for (instance_name, instance) in self._instances.items(): 425 237 426 # first, gather all unique variables referenced in any stage 427 # other than the last, independent of specific indices. this 428 # "gather" step is currently required because we're being lazy 429 # in terms of index management in the scenario tree - which 430 # should really be done in terms of sets of indices. 431 # NOTE: technically, the "instance variables" aren't really references 432 # to the variable in the instance - instead, the reference model. this 433 # isn't an issue now, but it could easily become one (esp. in avoiding deep copies). 434 instance_variables = {} 435 436 for stage in self._scenario_tree._stages[:-1]: 437 438 for (reference_variable, index_template, reference_indices) in stage._variables: 439 440 if reference_variable.name not in instance_variables.keys(): 441 442 instance_variables[reference_variable.name] = reference_variable 443 444 # for each blended variable, create a corresponding ph weight and average parameter in the instance. 445 # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor 446 # in the grand scheme of things. 447 448 for (variable_name, reference_variable) in instance_variables.items(): 449 450 # PH WEIGHT 451 452 new_w_index = reference_variable._index 453 new_w_parameter_name = "PHWEIGHT_"+reference_variable.name 454 # this bit of ugliness is due to Pyomo not correctly handling the Param construction 455 # case when the supplied index set consists strictly of None, i.e., the source variable 456 # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed. 457 new_w_parameter = None 458 if (len(new_w_index) is 1) and (None in new_w_index): 459 new_w_parameter = Param(name=new_w_parameter_name) 460 else: 461 new_w_parameter = Param(new_w_index,name=new_w_parameter_name) 462 setattr(instance,new_w_parameter_name,new_w_parameter) 463 464 # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference 465 # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not 466 # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the 467 # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo. 468 for index in new_w_index: 469 new_w_parameter[index] = 0.0 470 471 # PH AVG 472 473 new_avg_index = reference_variable._index 474 new_avg_parameter_name = "PHAVG_"+reference_variable.name 475 new_avg_parameter = None 476 if (len(new_avg_index) is 1) and (None in new_avg_index): 477 new_avg_parameter = Param(name=new_avg_parameter_name) 478 else: 479 new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name) 480 setattr(instance,new_avg_parameter_name,new_avg_parameter) 481 482 for index in new_avg_index: 483 new_avg_parameter[index] = 0.0 484 485 # PH RHO 486 487 new_rho_index = reference_variable._index 488 new_rho_parameter_name = "PHRHO_"+reference_variable.name 489 new_rho_parameter = None 490 if (len(new_avg_index) is 1) and (None in new_avg_index): 491 new_rho_parameter = Param(name=new_rho_parameter_name) 492 else: 493 new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name) 494 setattr(instance,new_rho_parameter_name,new_rho_parameter) 495 496 for index in new_rho_index: 497 new_rho_parameter[index] = self._rho 498 499 # PH LINEARIZED PENALTY TERM 500 501 if self._linearize_nonbinary_penalty_terms > 0: 502 503 new_penalty_term_variable_index = reference_variable._index 504 new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name 505 # this is a quadratic penalty term - the lower bound is 0! 506 new_penalty_term_variable = None 507 if (len(new_avg_index) is 1) and (None in new_avg_index): 508 new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None)) 509 else: 510 new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None)) 511 new_penalty_term_variable.construct() 512 setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable) 513 self._instance_augmented_attributes[instance_name].append(new_penalty_term_variable_name) 514 515 # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY. 516 517 # also controls whether weight updates proceed at any iteration. 518 519 new_blend_index = reference_variable._index 520 new_blend_parameter_name = "PHBLEND_"+reference_variable.name 521 new_blend_parameter = None 522 if (len(new_avg_index) is 1) and (None in new_avg_index): 523 new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary) 524 else: 525 new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary) 526 setattr(instance,new_blend_parameter_name,new_blend_parameter) 527 528 for index in new_blend_index: 529 new_blend_parameter[index] = 1 238 new_penalty_variable_names = create_ph_parameters(instance, self._scenario_tree, self._rho, self._linearize_nonbinary_penalty_terms) 239 self._instance_augmented_attributes[instance_name].extend(new_penalty_variable_names) 530 240 531 241 # a simple utility to extract the first-stage cost statistics, e.g., min, average, and max. … … 554 264 minimum_value = this_value 555 265 return minimum_value, sum_values/num_values, maximum_value 556 266 267 # a utility to transmit - across the PH solver manager - the current weights 268 # and averages for each of my problem instances. used to set up iteration K solves. 269 def _transmit_weights_and_averages(self): 270 271 for scenario_name, scenario_instance in self._instances.items(): 272 273 weights_to_transmit = [] 274 averages_to_transmit = [] 275 276 for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about. 277 278 for (variable, index_template, variable_indices) in stage._variables: 279 280 variable_name = variable.name 281 weight_parameter_name = "PHWEIGHT_"+variable_name 282 weights_to_transmit.append(getattr(scenario_instance, weight_parameter_name)) 283 average_parameter_name = "PHAVG_"+variable_name 284 averages_to_transmit.append(getattr(scenario_instance, average_parameter_name)) 285 286 self._solver_manager.transmit_weights_and_averages(scenario_instance, weights_to_transmit, averages_to_transmit) 287 288 # 289 # a utility to transmit - across the PH solver manager - the current rho values 290 # for each of my problem instances. mainly after PH iteration 0 is complete, 291 # in preparation for the iteration K solves. 292 # 293 294 def _transmit_rhos(self): 295 296 for scenario_name, scenario_instance in self._instances.items(): 297 298 rhos_to_transmit = [] 299 300 for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no rhos to worry about. 301 302 for (variable, index_template, variable_indices) in stage._variables: 303 304 variable_name = variable.name 305 rho_parameter_name = "PHRHO_"+variable_name 306 rhos_to_transmit.append(getattr(scenario_instance, rho_parameter_name)) 307 308 self._solver_manager.transmit_rhos(scenario_instance, rhos_to_transmit) 309 310 # 311 # a utility to transmit - across the PH solver manager - the current scenario 312 # tree node statistics to each of my problem instances. done prior to each 313 # PH iteration k. 314 # 315 316 def _transmit_tree_node_statistics(self): 317 318 for scenario_name, scenario_instance in self._instances.items(): 319 320 tree_node_minimums = {} 321 tree_node_maximums = {} 322 323 scenario = self._scenario_tree._scenario_map[scenario_name] 324 325 for tree_node in scenario._node_list: 326 327 tree_node_minimums[tree_node._name] = tree_node._minimums 328 tree_node_maximums[tree_node._name] = tree_node._maximums 329 330 self._solver_manager.transmit_tree_node_statistics(scenario_instance, tree_node_minimums, tree_node_maximums) 331 332 # 333 # a utility to enable - across the PH solver manager - weighted penalty objectives. 334 # 335 336 def _enable_ph_objectives(self): 337 338 for scenario_name, scenario_instance in self._instances.items(): 339 340 self._solver_manager.enable_ph_objective(scenario_instance) 341 557 342 558 343 """ Constructor … … 596 381 self._solver_manager_type = "serial" # serial or pyro are the options currently available 597 382 598 self._solver = None # will eventually be unnecessary once Bill eliminates the need for a solver in the solver manager constructor.383 self._solver = None 599 384 self._solver_manager = None 600 385 … … 1187 972 def form_iteration_k_objectives(self): 1188 973 1189 # for each blended variable (i.e., those not appearing in the final stage),1190 # add the linear and quadratic penalty terms to the objective.1191 974 for instance_name, instance in self._instances.items(): 1192 1193 objective_name = instance.active_components(Objective).keys()[0] 1194 objective = instance.active_components(Objective)[objective_name] 1195 # clone the objective, because we're going to augment (repeatedly) the original objective. 1196 objective_expression = self._original_objective_expression[instance_name].clone() 1197 # the quadratic expression is really treated as just a list - eventually should be treated as a full expression. 1198 quad_expression = 0.0 1199 1200 for stage in self._scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs 1201 1202 variable_tree_node = None 1203 for node in stage._tree_nodes: 1204 for scenario in node._scenarios: 1205 if scenario._name == instance_name: 1206 variable_tree_node = node 1207 break 1208 1209 for (reference_variable, index_template, variable_indices) in stage._variables: 1210 1211 variable_name = reference_variable.name 1212 variable_type = reference_variable.domain 1213 1214 w_parameter_name = "PHWEIGHT_"+variable_name 1215 w_parameter = instance.active_components(Param)[w_parameter_name] 1216 1217 average_parameter_name = "PHAVG_"+variable_name 1218 average_parameter = instance.active_components(Param)[average_parameter_name] 1219 1220 rho_parameter_name = "PHRHO_"+variable_name 1221 rho_parameter = instance.active_components(Param)[rho_parameter_name] 1222 1223 blend_parameter_name = "PHBLEND_"+variable_name 1224 blend_parameter = instance.active_components(Param)[blend_parameter_name] 1225 1226 node_min_parameter = variable_tree_node._minimums[variable_name] 1227 node_max_parameter = variable_tree_node._maximums[variable_name] 1228 1229 quad_penalty_term_variable = None 1230 if self._linearize_nonbinary_penalty_terms > 0: 1231 quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name 1232 quad_penalty_term_variable = instance.active_components(Var)[quad_penalty_term_variable_name] 1233 1234 instance_variable = instance.active_components(Var)[variable_name] 1235 1236 for index in variable_indices: 1237 1238 if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False): 1239 1240 # add the linear (w-weighted) term is a consistent fashion, independent of variable type. 1241 # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown. 1242 objective_expression += (w_parameter[index] * instance_variable[index]) 1243 1244 # there are some cases in which a user may want to avoid the proximal term completely. 1245 # it's of course only a good idea when there are at least bounds (both lb and ub) on 1246 # the variables to be blended. 1247 if self._drop_proximal_terms is False: 1248 1249 # deal with binaries 1250 if isinstance(variable_type, BooleanSet) is True: 1251 1252 if self._retain_quadratic_binary_terms is False: 1253 # this rather ugly form of the linearized quadratic expression term is required 1254 # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing 1255 # over the sum. 1256 new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \ 1257 (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \ 1258 (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index]) 1259 if objective.sense is minimize: 1260 objective_expression += new_term 1261 else: 1262 objective_expression -= new_term 1263 else: 1264 quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2) 1265 1266 # deal with everything else 1267 else: 1268 1269 if self._linearize_nonbinary_penalty_terms > 0: 1270 1271 # the variables are easy - just a single entry. 1272 if objective.sense is minimize: 1273 objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index]) 1274 else: 1275 objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index]) 1276 1277 # the constraints are somewhat nastier. 1278 1279 # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE" 1280 xavg = average_parameter[index] 1281 x = instance_variable[index] 1282 1283 lb = None 1284 ub = None 1285 1286 if x.lb is None: 1287 raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms" 1288 else: 1289 lb = x.lb() 1290 1291 if x.ub is None: 1292 raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms" 1293 else: 1294 ub = x.ub() 1295 1296 node_min = node_min_parameter[index]() 1297 node_max = node_max_parameter[index]() 1298 1299 # compute the breakpoint sequence according to the specified strategy. 1300 breakpoints = [] 1301 if self._breakpoint_strategy == 0: 1302 breakpoints = self.compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms) 1303 elif self._breakpoint_strategy == 1: 1304 breakpoints = self.compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms) 1305 elif self._breakpoint_strategy == 2: 1306 breakpoints = self.compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms) 1307 elif self._breakpoint_strategy == 3: 1308 breakpoints = self.compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms) 1309 else: 1310 raise ValueError, "A breakpoint distribution strategy="+str(self._breakpoint_strategy)+" is currently not supported within PH!" 1311 1312 for i in range(0,len(breakpoints)-1): 1313 1314 this_lb = breakpoints[i] 1315 this_ub = breakpoints[i+1] 1316 1317 piece_constraint_name = "QUAD_PENALTY_PIECE_"+str(i)+"_"+variable_name+str(index) 1318 if hasattr(instance, piece_constraint_name) is False: 1319 # this is the first time the constraint is being added - add it to the list of PH-specific constraints for this instance. 1320 self._instance_augmented_attributes[instance_name].append(piece_constraint_name) 1321 piece_constraint = Constraint(name=piece_constraint_name) 1322 piece_constraint.model = instance 1323 piece_expression = self._create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, quad_penalty_term_variable[index]) 1324 piece_constraint.add(None, piece_expression) 1325 setattr(instance, piece_constraint_name, piece_constraint) 1326 1327 else: 1328 1329 quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2) 1330 1331 # strictly speaking, this probably isn't necessary - parameter coefficients won't get 1332 # pre-processed out of the expression tree. however, if the under-the-hood should change, 1333 # we'll be covered. 1334 objective_expression.simplify(instance) 1335 instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression 1336 # if we are linearizing everything, then nothing will appear in the quadratic expression - 1337 # don't add the empty "0.0" expression to the objective. otherwise, the output file won't 1338 # be properly generated. 1339 if quad_expression != 0.0: 1340 instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression 975 976 new_attrs = form_ph_objective(instance_name, \ 977 instance, \ 978 self._original_objective_expression[instance_name], \ 979 self._scenario_tree, \ 980 self._linearize_nonbinary_penalty_terms, \ 981 self._drop_proximal_terms, \ 982 self._retain_quadratic_binary_terms, \ 983 self._breakpoint_strategy, \ 984 self._integer_tolerance) 985 self._instance_augmented_attributes[instance_name].extend(new_attrs) 1341 986 1342 987 def iteration_k_solve(self): … … 1350 995 1351 996 solve_start_time = time.time() 997 998 # STEP -1: if using a PH solver manager, propagate current weights/averages to the appropriate solver servers. 999 # ditto the tree node statistics, which are necessary if linearizing (so an optimization could be 1000 # performed here). 1001 # NOTE: We aren't currently propagating rhos, as they generally don't change - we need to 1002 # have a flag, though, indicating whether the rhos have changed, so they can be 1003 # transmitted if needed. 1004 if self._solver_manager_type == "ph": 1005 self._transmit_weights_and_averages() 1006 self._transmit_tree_node_statistics() 1352 1007 1353 1008 # STEP 0: set up all global solver options. … … 1443 1098 for objective_name in instance.active_components(Objective): 1444 1099 objective = instance.active_components(Objective)[objective_name] 1445 print "%20s %18.4f %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)1100 print "%20s %18.4f %14.4f" % (scenario._name, ph_objective_values[scenario._name], self._scenario_tree.compute_scenario_cost(instance)) 1446 1101 1447 1102 def solve(self): … … 1489 1144 for plugin in self._ph_plugins: 1490 1145 plugin.post_iteration_0(self) 1146 1147 # if using a PH solver server, trasnsmit the rhos prior to the iteration 1148 # k solve sequence. for now, we are assuming that the rhos don't change 1149 # on a per-iteration basis, but that assumption can be easily relaxed. 1150 # it is important to do this after the plugins have a chance to do their 1151 # computation. 1152 if self._solver_manager_type == "ph": 1153 self._transmit_rhos() 1154 self._enable_ph_objectives() 1491 1155 1492 1156 # checkpoint if it's time - which it always is after iteration 0, -
coopr.pysp/stable/2.3/coopr/pysp/phinit.py
r2361 r2434 531 531 print "Writing EF for remainder problem" 532 532 print "" 533 write_ef(ph._scenario_tree, ph._instances, os.path.expanduser(options.ef_output_file))533 create_and_write_ef(ph._scenario_tree, ph._instances, os.path.expanduser(options.ef_output_file)) 534 534 535 535 # … … 643 643 644 644 return ans 645 -
coopr.pysp/stable/2.3/coopr/pysp/phserver.py
r2361 r2434 6 6 # Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, 7 7 # the U.S. Government retains certain rights in this software. 8 # F or more information, see the Coopr README.txt file.8 # Fr more information, see the Coopr README.txt file. 9 9 # _________________________________________________________________________ 10 10 … … 23 23 from coopr.opt.base import SolverFactory 24 24 from coopr.pysp.scenariotree import * 25 from phutils import * 26 from phobjective import * 25 27 26 28 # Pyro 27 29 import Pyro.core 28 30 import Pyro.naming 31 from Pyro.errors import PyroError,NamingError 29 32 30 33 # garbage collection control. … … 35 38 import pstats 36 39 40 import pickle 41 42 # disable multi-threading. for a solver server, this 43 # would truly be a bad idea, as you (1) likely have limited 44 # solver licenses and (2) likely want to dedicate all compute 45 # resources on this node to a single solve. 46 Pyro.config.PYRO_MULTITHREADED=0 37 47 38 48 # … … 41 51 class PHSolverServer(object): 42 52 43 def __init__(self, scenario_instances, solver): 44 45 self._scenario_instances = scenario_instances 53 def __init__(self, scenario_instances, solver, scenario_tree): 54 55 # the obvious stuff! 56 self._instances = scenario_instances # a map from scenario name to pyomo model instance 46 57 self._solver = solver 58 self._scenario_tree = scenario_tree 59 60 # the objective functions are modified throughout the course of PH iterations. 61 # save the original, as a baseline to modify in subsequent iterations. reserve 62 # the original objectives, for subsequent modification. 63 self._original_objective_expression = {} 64 for instance_name, instance in self._instances.items(): 65 objective_name = instance.active_components(Objective).keys()[0] 66 expr = instance.active_components(Objective)[objective_name]._data[None].expr 67 if isinstance(expr, Expression) is False: 68 expr = _IdentityExpression(expr) 69 self._original_objective_expression[instance_name] = expr 70 71 # the server is in one of two modes - solve the baseline instances, or 72 # those augmented with the PH penalty terms. default is standard. 73 # NOTE: This is currently a global flag for all scenarios handled 74 # by this server - easy enough to extend if we want. 75 self._solving_standard_objective = True 76 77 def enable_standard_objective(self): 78 79 # print "RECEIVED REQUEST TO ENABLE STANDARD OBJECTIVE FOR ALL SCENARIOS" 80 81 self._solving_standard_objective = True 82 83 def enable_ph_objective(self): 84 85 # print "RECEIVED REQUEST TO ENABLE PH OBJECTIVE FOR ALL SCENARIOS" 86 87 self._solving_standard_objective = False 47 88 48 89 def solve(self, scenario_name): 49 90 50 print "RECEIVED SOLVE REQUEST!" 51 print "SCENARIO TO SOLVE=",scenario_name 52 if scenario_name not in self._scenario_instances: 91 # print "RECEIVED REQUEST TO SOLVE SCENARIO INSTANCE="+scenario_name 92 # if self._solving_standard_objective is True: 93 # print "OBJECTIVE=STANDARD" 94 # else: 95 # print "OBJECTIVE=PH" 96 97 if scenario_name not in self._instances: 53 98 print "***ERROR: Requested instance to solve not in PH solver server instance collection!" 54 99 return None 55 scenario_instance = self._scenario_instances[scenario_name] 56 print "SCENARIO INSTANCE=",scenario_instance 57 print "SUPPORTS MULTITHREADING=",Pyro.util.supports_multithreading() 58 print "SOLVING" 59 self._solver.solve(scenario_instance) 60 print "DONE WITH SOLVE" 61 return scenario_instance 62 # return result 63 100 scenario_instance = self._instances[scenario_name] 101 102 # form the desired objective, depending on the solve mode. 103 if self._solving_standard_objective is True: 104 form_standard_objective(scenario_name, scenario_instance, self._original_objective_expression[scenario_name], self._scenario_tree) 105 else: 106 # TBD - command-line drive the various options (integer tolerance, breakpoint strategies, etc.) 107 form_ph_objective(scenario_name, scenario_instance, \ 108 self._original_objective_expression[scenario_name], self._scenario_tree, \ 109 False, False, False, 0, 0.00001) 110 111 # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you 112 # don't do this, you won't see any chance in the output files as you vary the problem parameters! 113 # ditto for instance fixing! 114 scenario_instance.preprocess() 115 116 results = self._solver.solve(scenario_instance) 117 118 print "Successfully solved scenario instance="+scenario_name 119 # print "RESULTS:" 120 # results.write(num=1) 121 encoded_results = pickle.dumps(results) 122 123 return encoded_results 124 125 def update_weights_and_averages(self, scenario_name, new_weights, new_averages): 126 127 # print "RECEIVED REQUEST TO UPDATE WEIGHTS AND AVERAGES FOR SCENARIO=",scenario_name 128 129 if scenario_name not in self._instances: 130 print "***ERROR: Received request to update weights for instance not in PH solver server instance collection!" 131 return None 132 scenario_instance = self._instances[scenario_name] 133 134 for weight_update in new_weights: 135 136 weight_index = weight_update._index 137 138 target_weight_parameter = getattr(scenario_instance, weight_update.name) 139 140 for index in weight_index: 141 target_weight_parameter[index] = weight_update[index]() 142 143 for average_update in new_averages: 144 145 average_index = average_update._index 146 147 target_average_parameter = getattr(scenario_instance, average_update.name) 148 149 for index in average_index: 150 target_average_parameter[index] = average_update[index]() 151 152 def update_rhos(self, scenario_name, new_rhos): 153 154 # print "RECEIVED REQUEST TO UPDATE RHOS FOR SCENARIO=",scenario_name 155 156 if scenario_name not in self._instances: 157 print "***ERROR: Received request to update weights for instance not in PH solver server instance collection!" 158 return None 159 scenario_instance = self._instances[scenario_name] 160 161 for rho_update in new_rhos: 162 163 rho_index = rho_update._index 164 165 target_rho_parameter = getattr(scenario_instance, rho_update.name) 166 167 for index in rho_index: 168 target_rho_parameter[index] = rho_update[index]() # the value operator is crucial! 169 170 def update_tree_node_statistics(self, scenario_name, new_node_minimums, new_node_maximums): 171 172 # print "RECEIVED REQUEST TO UPDATE TREE NODE STATISTICS SCENARIO=",scenario_name 173 174 for tree_node_name, tree_node_minimums in new_node_minimums.items(): 175 176 tree_node = self._scenario_tree._tree_node_map[tree_node_name] 177 tree_node._minimums = tree_node_minimums 178 179 for tree_node_name, tree_node_maximums in new_node_maximums.items(): 180 181 tree_node = self._scenario_tree._tree_node_map[tree_node_name] 182 tree_node._maximums = tree_node_maximums 183 184 64 185 # 65 186 # utility method to construct an option parser for ph arguments, to be … … 80 201 dest="scenarios", 81 202 default=[]) 203 parser.add_option("--all-scenarios", 204 help="Indicate that the server is responsible for solving all scenarios", 205 action="store_true", 206 dest="all_scenarios", 207 default=False) 82 208 parser.add_option("--report-solutions", 83 209 help="Always report PH solutions after each iteration. Enabled if --verbose is enabled. Default is False.", … … 151 277 type="int", 152 278 default=0) 153 parser.add_option("-- enable-gc",154 help=" Enable the python garbage collecter. Default is True.",155 action="store_true", 156 dest=" enable_gc",157 default= True)279 parser.add_option("--disable-gc", 280 help="Disable the python garbage collecter. Default is False.", 281 action="store_true", 282 dest="disable_gc", 283 default=False) 158 284 159 285 parser.usage=usage_string … … 182 308 print "" 183 309 print "***ERROR: Exiting test driver: No 'model' object created in module "+reference_model_filename 184 return None, None, None 310 return None, None, None, None 185 311 186 312 if model_import.model is None: 187 313 print "" 188 314 print "***ERROR: Exiting test driver: 'model' object equals 'None' in module "+reference_model_filename 189 return None, None, None 315 return None, None, None, None 190 316 191 317 reference_model = model_import.model 192 318 except IOError: 193 319 print "***ERROR: Failed to load scenario reference model from file="+reference_model_filename 194 return None, None, None 320 return None, None, None, None 195 321 196 322 try: … … 201 327 except IOError: 202 328 print "***ERROR: Failed to load scenario reference instance data from file="+reference_instance_filename 203 return None, None, None 329 return None, None, None, None 204 330 205 331 # … … 231 357 # 232 358 233 def run_server(options, scenario_instances, solver ):359 def run_server(options, scenario_instances, solver, scenario_tree): 234 360 235 361 start_time = time.time() 236 362 237 print "RUNNING SERVER NOW!" 238 Pyro.core.initServer(banner=1) 239 240 print "NUMBER OF SCENARIO INSTANCES="+str(len(scenario_instances)) 241 242 print "ATTEMPTING TO LOCATE NAME SERVER" 363 Pyro.core.initServer(banner=0) 364 243 365 locator = Pyro.naming.NameServerLocator() 244 ns = locator.getNS() 245 print "FOUND NAME SERVER" 366 ns = None 367 try: 368 ns = locator.getNS() 369 except Pyro.errors.NamingError: 370 raise RuntimeError, "PH solver server failed to locate Pyro name server" 246 371 247 372 solver_daemon = Pyro.core.Daemon() 248 373 solver_daemon.useNameServer(ns) 249 374 250 daemon_object = PHSolverServer(scenario_instances, solver )375 daemon_object = PHSolverServer(scenario_instances, solver, scenario_tree) 251 376 delegator_object = Pyro.core.ObjBase() 252 377 delegator_object.delegateTo(daemon_object) … … 255 380 # each scenario processed by this daemon. 256 381 for scenario_name, scenario_instance in scenario_instances.items(): 257 solver_daemon.connect(delegator_object, scenario_name) 258 259 while True: 260 solver_daemon.handleRequests() 261 print "HANDLED A REQUEST!" 262 263 solver_daemon.disconnect() 264 265 print "DONE RUNNING SERVER!" 382 # first, see if it's already registered - if no, cull the old entry. 383 # NOTE: This is a hack, as the cleanup / disconnect code doesn't 384 # seem to be invoked when the daemon is killed. 385 try: 386 ns.unregister(scenario_name) 387 except NamingError: 388 pass 389 390 try: 391 solver_daemon.connect(delegator_object, scenario_name) 392 except Pyro.errors.NamingError: 393 raise RuntimeError, "Entry in name server already exists for scenario="+scenario_name 394 395 try: 396 solver_daemon.requestLoop() 397 except KeyboardInterrupt, SystemExit: 398 pass 399 400 solver_daemon.disconnect(delegator_object) 401 solver_daemon.shutdown() 266 402 267 403 end_time = time.time() … … 273 409 def exec_server(options): 274 410 275 # we need the base model (not so much the reference instance, but that -276 # can't hurt too much until proven otherwise, that is) to construct411 # we need the base model (not so much the reference instance, but that 412 # can't hurt too much - until proven otherwise, that is) to construct 277 413 # the scenarios that this server is responsible for. 278 print "LOADING REFERENCE AND SCENARIO MODELS/INSTANCES"279 414 reference_model, reference_instance, scenario_tree, scenario_tree_instance = load_reference_and_scenario_models(options) 280 415 if reference_model is None or reference_instance is None or scenario_tree is None: … … 284 419 # construct the set of scenario instances based on the command-line. 285 420 scenario_instances = {} 286 if len(options.scenarios) == 0: 287 print "***Unable to launch PH solver server - no scenario specified!" 421 if (options.all_scenarios is True) and (len(options.scenarios) > 0): 422 print "***WARNING: Both individual scenarios and all-scenarios were specified on the command-line; proceeding using all scenarios" 423 424 if (len(options.scenarios) == 0) and (options.all_scenarios is False): 425 print "***Unable to launch PH solver server - no scenario(s) specified!" 288 426 return 289 427 290 for scenario_name in options.scenarios: 428 scenarios_to_construct = [] 429 if options.all_scenarios is True: 430 scenarios_to_construct.extend(scenario_tree._scenario_map.keys()) 431 else: 432 scenarios_to_construct.extend(options.scenarios) 433 434 for scenario_name in scenarios_to_construct: 291 435 print "Creating instance for scenario="+scenario_name 292 436 if scenario_tree.contains_scenario(scenario_name) is False: … … 294 438 return 295 439 440 # create the baseline scenario instance 296 441 scenario_instance = construct_scenario_instance(scenario_tree, 297 442 options.instance_directory, … … 301 446 if scenario_instance is None: 302 447 print "***Unable to launch PH solver server - failed to create instance for scenario="+scenario_name 448 return 449 450 # augment the instance with PH-specific parameters (weights, rhos, etc). 451 # TBD: The default rho of 1.0 is kind of bogus. Need to somehow propagate 452 # this value and the linearization parameter as a command-line argument. 453 new_penalty_variable_names = create_ph_parameters(scenario_instance, scenario_tree, 1.0, False) 454 303 455 scenario_instances[scenario_name] = scenario_instance 304 456 … … 318 470 solver._report_timing = True 319 471 320 print "SOLVING!"321 results = solver.solve(scenario_instance)322 print "TYPE OF RESULTS=",type(results)323 results.write(num=1)324 # print "RESULTS=",results325 print "DONE SOLVING!"326 327 472 # spawn the daemon. 328 run_server(options, scenario_instances, solver )473 run_server(options, scenario_instances, solver, scenario_tree) 329 474 330 475 def run(args=None): … … 349 494 # much sense - so it is disabled by default. Because: It drops 350 495 # the run-time by a factor of 3-4 on bigger instances. 351 if options. enable_gc is False:496 if options.disable_gc is True: 352 497 gc.disable() 353 498 else: -
coopr.pysp/stable/2.3/coopr/pysp/phutils.py
r2391 r2434 224 224 225 225 return scenario_instance 226 227 228 # creates all PH parameters for a problem instance, given a scenario tree 229 # (to identify the relevant variables), a default rho (simply for initialization), 230 # and a boolean indicating if quadratic penalty terms are to be linearized. 231 # returns a list of any created variables, specifically when linearizing - 232 # this is useful to clean-up reporting. 233 234 def create_ph_parameters(instance, scenario_tree, default_rho, linearizing_penalty_terms): 235 236 new_penalty_variable_names = [] 237 238 # first, gather all unique variables referenced in any stage 239 # other than the last, independent of specific indices. this 240 # "gather" step is currently required because we're being lazy 241 # in terms of index management in the scenario tree - which 242 # should really be done in terms of sets of indices. 243 # NOTE: technically, the "instance variables" aren't really references 244 # to the variable in the instance - instead, the reference model. this 245 # isn't an issue now, but it could easily become one (esp. in avoiding deep copies). 246 instance_variables = {} 247 248 for stage in scenario_tree._stages[:-1]: 249 250 for (reference_variable, index_template, reference_indices) in stage._variables: 251 252 if reference_variable.name not in instance_variables.keys(): 253 254 instance_variables[reference_variable.name] = reference_variable 255 256 # for each blended variable, create a corresponding ph weight and average parameter in the instance. 257 # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor 258 # in the grand scheme of things. 259 260 for (variable_name, reference_variable) in instance_variables.items(): 261 262 # PH WEIGHT 263 264 new_w_index = reference_variable._index 265 new_w_parameter_name = "PHWEIGHT_"+reference_variable.name 266 # this bit of ugliness is due to Pyomo not correctly handling the Param construction 267 # case when the supplied index set consists strictly of None, i.e., the source variable 268 # is a singleton. this case be cleaned up when the source issue in Pyomo is fixed. 269 new_w_parameter = None 270 if (len(new_w_index) is 1) and (None in new_w_index): 271 new_w_parameter = Param(name=new_w_parameter_name) 272 else: 273 new_w_parameter = Param(new_w_index,name=new_w_parameter_name) 274 setattr(instance,new_w_parameter_name,new_w_parameter) 275 276 # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference 277 # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not 278 # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the 279 # issue in any case for now. NOTE: I'm not sure this is still the case in Pyomo. 280 for index in new_w_index: 281 new_w_parameter[index] = 0.0 282 283 # PH AVG 284 285 new_avg_index = reference_variable._index 286 new_avg_parameter_name = "PHAVG_"+reference_variable.name 287 new_avg_parameter = None 288 if (len(new_avg_index) is 1) and (None in new_avg_index): 289 new_avg_parameter = Param(name=new_avg_parameter_name) 290 else: 291 new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name) 292 setattr(instance,new_avg_parameter_name,new_avg_parameter) 293 294 for index in new_avg_index: 295 new_avg_parameter[index] = 0.0 296 297 # PH RHO 298 299 new_rho_index = reference_variable._index 300 new_rho_parameter_name = "PHRHO_"+reference_variable.name 301 new_rho_parameter = None 302 if (len(new_avg_index) is 1) and (None in new_avg_index): 303 new_rho_parameter = Param(name=new_rho_parameter_name) 304 else: 305 new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name) 306 setattr(instance,new_rho_parameter_name,new_rho_parameter) 307 308 for index in new_rho_index: 309 new_rho_parameter[index] = default_rho 310 311 # PH LINEARIZED PENALTY TERM 312 313 if linearizing_penalty_terms > 0: 314 315 new_penalty_term_variable_index = reference_variable._index 316 new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name 317 # this is a quadratic penalty term - the lower bound is 0! 318 new_penalty_term_variable = None 319 if (len(new_avg_index) is 1) and (None in new_avg_index): 320 new_penalty_term_variable = Var(name=new_penalty_term_variable_name, bounds=(0.0,None)) 321 else: 322 new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None)) 323 new_penalty_term_variable.construct() 324 setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable) 325 new_penalty_variable_names.append(new_penalty_term_variable_name) 326 327 # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY. 328 329 # also controls whether weight updates proceed at any iteration. 330 331 new_blend_index = reference_variable._index 332 new_blend_parameter_name = "PHBLEND_"+reference_variable.name 333 new_blend_parameter = None 334 if (len(new_avg_index) is 1) and (None in new_avg_index): 335 new_blend_parameter = Param(name=new_blend_parameter_name, within=Binary) 336 else: 337 new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary) 338 setattr(instance,new_blend_parameter_name,new_blend_parameter) 339 340 for index in new_blend_index: 341 new_blend_parameter[index] = 1 342 343 return new_penalty_variable_names -
coopr.pysp/stable/2.3/coopr/pysp/scenariotree.py
r2391 r2434 11 11 import sys 12 12 import types 13 from coopr.pyomo import *14 13 import copy 15 14 import os.path 16 15 import traceback 17 16 17 from coopr.pyomo import * 18 18 from phutils import * 19 19 … … 23 23 24 24 """ 25 def __init__(self, *args, **kwds):26 27 self._name = ""28 self._stage = None25 def __init__(self, name, conditional_probability, stage, reference_instance): 26 27 self._name = name 28 self._stage = stage 29 29 self._parent = None 30 30 self._children = [] # a collection of ScenarioTreeNodes 31 self._conditional_probability = None# conditional on parent31 self._conditional_probability = conditional_probability # conditional on parent 32 32 self._scenarios = [] # a collection of all Scenarios passing through this node in the tree 33 33 … … 44 44 self._maximums = {} 45 45 46 # solution (variable) values for this node. assumed to be distinct 47 # from self._averages, as the latter are not necessarily feasible. 48 # objects in the map are actual variables. 49 self._solution = {} 50 51 # for each variable referenced in the stage, clone the variable 52 # for purposes of storing solutions. we are being wasteful in 53 # terms copying indices that may not be referenced in the stage. 54 # this is something that we might revisit if space/performance 55 # is an issue (space is the most likely issue) 56 for variable, match_template, variable_indices in self._stage._variables: 57 58 # don't bother copying bounds for variables, as the values stored 59 # here are computed elsewhere - and that entity is responsible for 60 # ensuring feasibility. this also leaves room for specifying infeasible 61 # or partial solutions. 62 63 new_variable_index = variable._index 64 new_variable_name = variable._name 65 new_variable = None 66 if (len(new_variable_index) is 1) and (None in new_variable_index): 67 new_variable = Var(name=new_variable_name) 68 else: 69 new_variable = Var(new_variable_index, name=new_variable_name) 70 71 self._solution[new_variable_name] = new_variable 72 46 73 # 47 74 # a utility to compute the cost of the current node plus the expected costs of child nodes. 48 75 # 76 49 77 def computeExpectedNodeCost(self, scenario_instance_map): 50 78 … … 90 118 """ Constructor 91 119 Arguments: 92 scenarioinstance - the referencescenario instance.120 scenarioinstance - the reference (deterministic) scenario instance. 93 121 scenariotreeinstance - the pyomo model specifying all scenario tree (text) data. 94 122 scenariobundlelist - a list of scenario names to retain, i.e., cull the rest to create a reduced tree! … … 172 200 # can't do a single pass because the objects may not exist. 173 201 for tree_node_name in node_ids: 174 new_tree_node = ScenarioTreeNode() 175 new_tree_node._name = tree_node_name 202 203 if tree_node_name not in node_stage_ids: 204 raise ValueError, "No stage is assigned to tree node=" + tree_node._name 205 206 stage_name = node_stage_ids[tree_node_name].value 207 if stage_name not in self._stage_map.keys(): 208 raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name 209 210 new_tree_node = ScenarioTreeNode(tree_node_name, 211 node_probability_map[tree_node_name].value, 212 self._stage_map[stage_name], 213 self._reference_instance) 176 214 177 215 self._tree_nodes.append(new_tree_node) 178 216 self._tree_node_map[tree_node_name] = new_tree_node 179 180 new_tree_node._conditional_probability = node_probability_map[tree_node_name].value 181 182 if tree_node_name not in node_stage_ids: 183 raise ValueError, "No stage is assigned to tree node=" + tree_node._name 184 else: 185 stage_name = node_stage_ids[new_tree_node._name].value 186 if stage_name not in self._stage_map.keys(): 187 raise ValueError, "Unknown stage=" + stage_name + " assigned to tree node=" + tree_node._name 188 else: 189 new_tree_node._stage = self._stage_map[stage_name] 190 self._stage_map[stage_name]._tree_nodes.append(new_tree_node) 217 self._stage_map[stage_name]._tree_nodes.append(new_tree_node) 191 218 192 219 # link up the tree nodes objects based on the child id sets. … … 347 374 # is the indicated scenario in the tree? 348 375 # 376 349 377 def contains_scenario(self, name): 350 378 return name in self._scenario_map.keys() … … 353 381 # get the scenario object from the tree. 354 382 # 383 355 384 def get_scenario(self, name): 356 385 return self._scenario_map[name] 386 387 # 388 # compute the scenario cost for the input instance, i.e., 389 # the sum of all stage cost variables. 390 # 391 392 def compute_scenario_cost(self, instance): 393 aggregate_cost = 0.0 394 for stage in self._stages: 395 instance_cost_variable = instance.active_components(Var)[stage._cost_variable[0].name][stage._cost_variable[1]]() 396 aggregate_cost += instance_cost_variable 397 return aggregate_cost 357 398 358 399 # … … 363 404 # of the scenario tree structure. 364 405 # 406 365 407 def compress(self, scenario_bundle_list): 366 408 … … 450 492 # returns the root node of the scenario tree 451 493 # 494 452 495 def findRootNode(self): 453 496 … … 461 504 # the maximal length of identifiers in various categories. 462 505 # 506 463 507 def computeIdentifierMaxLengths(self): 464 508 … … 471 515 # a utility function to (partially, at the moment) validate a scenario tree 472 516 # 517 473 518 def validate(self): 474 519 -
coopr.pysp/stable/2.3/coopr/pysp/tests/unit/test_ph.py
r2391 r2434 18 18 # Import the testing packages 19 19 # 20 import unittest21 20 import pyutilib.misc 22 import pyutilib.th 21 import pyutilib.th as unittest 23 22 import pyutilib.subprocess 24 23 import coopr.pysp … … 46 45 47 46 # 48 # Define a testing class, using the pyutilib.th.TestCase class. This is 49 # an extension of the unittest.TestCase class that adds additional testing 50 # functions. 47 # Define a testing class, using the unittest.TestCase class. 51 48 # 52 class TestPH( pyutilib.th.TestCase):49 class TestPH(unittest.TestCase): 53 50 54 51 def cleanup(self): … … 61 58 del sys.modules["ReferenceModel"] 62 59 60 @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available") 63 61 def test_quadratic_farmer(self): 64 65 if cplex_available is False:66 return67 68 62 pyutilib.misc.setup_redirect(current_directory+"farmer_quadratic.out") 69 63 farmer_examples_dir = pysp_examples_dir + "farmer" … … 77 71 self.failUnlessFileEqualsBaseline(current_directory+"farmer_quadratic.out",current_directory+"farmer_quadratic.baseline", filter=filter_time) 78 72 73 @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available") 79 74 def test_linearized_farmer(self): 80 81 solver_string = ""82 75 if cplex_available: 83 76 solver_string="cplex" 84 #elif glpk_available:85 #solver_string="glpk"86 else:87 return88 89 77 pyutilib.misc.setup_redirect(current_directory+"farmer_linearized.out") 90 78 farmer_examples_dir = pysp_examples_dir + "farmer" … … 98 86 self.failUnlessFileEqualsBaseline(current_directory+"farmer_linearized.out",current_directory+"farmer_linearized.baseline", filter=filter_time) 99 87 88 @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available") 100 89 def test_linearized_farmer_nodedata(self): 101 102 solver_string = ""103 90 if cplex_available: 104 91 solver_string="cplex" 105 #elif glpk_available:106 #solver_string="glpk"107 else:108 return109 110 92 pyutilib.misc.setup_redirect(current_directory+"farmer_linearized_nodedata.out") 111 93 farmer_examples_dir = pysp_examples_dir + "farmer" … … 119 101 self.failUnlessFileEqualsBaseline(current_directory+"farmer_linearized_nodedata.out",current_directory+"farmer_linearized_nodedata.baseline", filter=filter_time) 120 102 103 @unittest.skipIf(not cplex_available, "The 'cplex' executable is not available") 121 104 def test_quadratic_sizes3(self): 122 105 123 106 # Ignore this for now 124 107 return 125 126 if cplex_available is False:127 return128 108 129 109 pyutilib.misc.setup_redirect(current_directory+"sizes3_quadratic.out") … … 146 126 147 127 if __name__ == "__main__": 148 unittest.main()128 unittest.main() -
coopr.pysp/stable/2.3/scripts/ph_test_client
r2361 r2434 4 4 import Pyro.naming 5 5 6 from coopr.opt.base import SolverFactory 7 from coopr.opt import * 6 import pickle 7 8 import pyutilib.misc 9 import pyutilib.component.core 10 from coopr.opt.results import SolverResults 8 11 9 12 import sys … … 38 41 39 42 print "CALLING SOLVE METHOD" 40 result = obj.solve(scenario_name) 43 encoded_result = obj.solve(scenario_name) 44 result = pickle.loads(encoded_result) 45 print "SUCCESSFULLY OBTAINED RESULT!" 41 46 42 print "RESULT=",result 47 print "RESULTS:" 48 result.write(num=1) 43 49 44 50 print "DONE"
Note: See TracChangeset
for help on using the changeset viewer.