source: coopr.pysp/trunk/coopr/pysp/ef.py @ 2110

Last change on this file since 2110 was 2110, checked in by jwatson, 10 years ago

Promoted the scenario tree model to a pysp/util package (with Bill's excellent guidance), and eliminated the replicated ScenarioStructure?.py files from the examples directory.

The only place the master scenario tree model exists is now in coopr/pysp/util/scenariomodels.py, and the (abstract) model is called scenario_tree_model.

File size: 41.8 KB
Line 
1import pyutilib
2import sys
3import os
4import time
5import traceback
6import copy
7
8from coopr.pysp.scenariotree import *
9from coopr.pysp.convergence import *
10from coopr.pysp.ph import *
11
12from coopr.pyomo.base import *
13from coopr.pyomo.io import *
14
15from coopr.pyomo.base.var import _VarValue, _VarBase
16
17# brain-dead utility for determing if there is a binary to write in the
18# composite model - need to know this, because CPLEX doesn't like empty
19# binary blocks in the LP file.
20def binaries_present(master_model, scenario_instances):
21
22   # check the master model first.
23   for var in master_model.active_components(Var).values():
24      if isinstance(var.domain, BooleanSet):
25         return True
26
27   # scan the scenario instances next.
28   for scenario_name in scenario_instances.keys():
29      scenario_instance = scenario_instances[scenario_name]
30      for var in scenario_instance.active_components(Var).values():
31         if isinstance(var.domain, BooleanSet):
32            return True
33
34   return False
35
36# brain-dead utility for determing if there is a binary to write in the
37# composite model - need to know this, because CPLEX doesn't like empty
38# integer blocks in the LP file.
39def integers_present(master_model, scenario_instances):
40
41   # check the master model first.
42   for var in master_model.active_components(Var).values():
43      if isinstance(var.domain, IntegerSet):
44         return True
45
46   # scan the scenario instances next.
47   for scenario_name in scenario_instances.keys():
48      scenario_instance = scenario_instances[scenario_name]
49      for var in scenario_instance.active_components(Var).values():
50         if isinstance(var.domain, IntegerSet):
51            return True
52
53   return False
54
55# the main extensive-form writer routine - including read of scenarios/etc.
56def write_ef_from_scratch(model_directory, instance_directory, output_filename, \
57                          generate_weighted_cvar, cvar_weight, risk_alpha):
58
59   start_time = time.time()
60
61   scenario_data_directory_name = instance_directory
62
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)
82      print ""
83   
84   #
85   # create and populate the core model
86   #
87   master_scenario_model = None
88   master_scenario_instance = None
89   
90   try:
91     
92      reference_model_filename = model_directory+os.sep+"ReferenceModel.py"   
93      modelimport = pyutilib.misc.import_file(reference_model_filename)
94      if "model" not in dir(modelimport):
95         print ""
96         print "Exiting ef module: No 'model' object created in module "+reference_model_filename
97         sys.exit(0)
98      if modelimport.model is None:
99         print ""
100         print "Exiting ef module: 'model' object equals 'None' in module "+reference_model_filename
101         sys.exit(0)
102   
103      master_scenario_model = modelimport.model
104
105   except IOError:
106     
107      print "***ERROR: Failed to load scenario reference model from file="+reference_model_filename
108      return
109
110   try:
111     
112      reference_scenario_filename = instance_directory+os.sep+"ReferenceModel.dat"
113      master_scenario_instance = master_scenario_model.create(reference_scenario_filename)
114     
115   except IOError:
116     
117      print "***ERROR: Failed to load scenario reference instance data from file="+reference_instance_filename
118      return           
119
120   #
121   # create and populate the scenario tree model
122   #
123
124   from coopr.pysp.util.scenariomodels import scenario_tree_model
125
126   tree_data = scenario_tree_model.create(instance_directory+os.sep+"ScenarioStructure.dat")
127
128   #
129   # construct the scenario tree
130   #
131   scenario_tree = ScenarioTree(model=master_scenario_instance,
132                                nodes=tree_data.Nodes,
133                                nodechildren=tree_data.Children,
134                                nodestages=tree_data.NodeStage,
135                                nodeprobabilities=tree_data.ConditionalProbability,
136                                stages=tree_data.Stages,
137                                stagevariables=tree_data.StageVariables,
138                                stagecostvariables=tree_data.StageCostVariable,
139                                scenarios=tree_data.Scenarios,
140                                scenarioleafs=tree_data.ScenarioLeafNode,
141                                scenariobaseddata=tree_data.ScenarioBasedData)
142
143   #
144   # print the input tree for validation/information purposes.
145   #
146   scenario_tree.pprint()
147
148   #
149   # validate the tree prior to doing anything serious
150   #
151   print ""
152   if scenario_tree.validate() is False:
153      print "***Scenario tree is invalid****"
154      sys.exit(1)
155   else:
156      print "Scenario tree is valid!"
157   print ""
158
159   #
160   # construct instances for each scenario
161   #
162
163   instances = {}
164   
165   if scenario_tree._scenario_based_data == 1:
166      print "Scenario-based instance initialization enabled"
167   else:
168      print "Node-based instance initialization enabled"
169         
170   for scenario in scenario_tree._scenarios:
171
172      scenario_instance = None
173
174      print "Creating instance for scenario=" + scenario._name
175
176      try:
177         if scenario_tree._scenario_based_data == 1:
178            scenario_data_filename = scenario_data_directory_name + os.sep + scenario._name + ".dat"
179#            print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
180            scenario_instance = master_scenario_model.create(scenario_data_filename)
181         else:
182            scenario_instance = master_scenario_model.clone()
183            scenario_data = ModelData()
184            current_node = scenario._leaf_node
185            while current_node is not None:
186               node_data_filename = scenario_data_directory_name + os.sep + current_node._name + ".dat"
187#               print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
188               scenario_data.add_data_file(node_data_filename)
189               current_node = current_node._parent
190            scenario_data.read(model=scenario_instance)
191            scenario_instance._load_model_data(scenario_data)
192            scenario_instance.presolve()
193      except:
194         print "Encountered exception in model instance creation - traceback:"
195         traceback.print_exc()
196         raise RuntimeError, "Failed to create model instance for scenario=" + scenario._name
197
198      # name each instance with the scenario name, so the prefixes in the EF make sense.
199      scenario_instance.name = scenario._name
200     
201      scenario_instance.presolve()
202      instances[scenario._name] = scenario_instance
203
204   print ""
205
206   ################################################################################################
207   #### CREATE THE MASTER / BINDING INSTANCE ######################################################
208   ################################################################################################
209
210   master_binding_instance = Model()
211   master_binding_instance.name = "MASTER"
212
213   # walk the scenario tree - create variables representing the common values for all scenarios
214   # associated with that node. the constraints will be created later. also create expected-cost
215   # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
216   # master variables are created for all nodes but those in the last stage. expected cost variables
217   # are, for no particularly good reason other than easy coding, created for nodes in all stages.
218   print "Creating variables for master binding instance"
219
220   for stage in scenario_tree._stages:
221
222      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
223
224         print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=", stage_variable_indices
225
226         for tree_node in stage._tree_nodes:
227
228            if stage != scenario_tree._stages[-1]:     
229
230               master_variable_name = tree_node._name + "_" + stage_variable.name
231
232               # because there may be a single stage variable and multiple indices, check
233               # for the existence of the variable at this node - if you don't, you'll
234               # inadvertently over-write what was there previously!
235               master_variable = None
236               try:
237                  master_variable = getattr(master_binding_instance, master_variable_name)
238               except:
239                  # the deepcopy is probably too expensive (and unnecessary) computationally -
240                  # easier to just use the constructor with the stage variable index/bounds/etc.
241                  # NOTE: need to re-assign the master variables for each _varval - they probably
242                  #       point to a bogus model.
243                  new_master_variable = copy.deepcopy(stage_variable)
244                  new_master_variable.name = master_variable_name
245                  new_master_variable._model = master_binding_instance
246                  setattr(master_binding_instance, master_variable_name, new_master_variable)
247
248                  master_variable = new_master_variable
249
250               for index in stage_variable_indices:
251
252                  is_used = True # until proven otherwise                     
253                  for scenario in tree_node._scenarios:
254                     instance = instances[scenario._name]
255                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
256                        is_used = False
257
258                  is_fixed = False # until proven otherwise
259                  for scenario in tree_node._scenarios:
260                     instance = instances[scenario._name]
261                     if getattr(instance,stage_variable.name)[index].fixed is True:
262                        is_fixed = True
263
264                  if (is_used is True) and (is_fixed is False):
265                           
266                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
267                     # and because presolve/simplification is name-based, the names *have* to be different.
268                     master_variable[index].var = master_variable
269                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
270
271                     for scenario in tree_node._scenarios:
272
273                        scenario_instance = instances[scenario._name]
274                        scenario_variable = getattr(scenario_instance, stage_variable.name)
275                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
276                        new_constraint = Constraint(name=new_constraint_name)
277                        new_expr = master_variable[index] - scenario_variable[index]
278                        new_constraint.add(None, (0.0, new_expr, 0.0))
279                        new_constraint._model = master_binding_instance
280                        setattr(master_binding_instance, new_constraint_name, new_constraint)
281
282            # create a variable to represent the expected cost at this node -
283            # the constraint to compute this comes later.
284            expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
285            expected_cost_variable = Var(name=expected_cost_variable_name)
286            expected_cost_variable._model = master_binding_instance
287            setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
288
289   # if we're generating the weighted CVaR objective term, create the corresponding variable and
290   # the master CVaR eta variable.
291   if generate_weighted_cvar is True:
292      root_node = scenario_tree._stages[0]._tree_nodes[0]
293     
294      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
295      cvar_cost_variable = Var(name=cvar_cost_variable_name)
296      setattr(master_binding_instance, cvar_cost_variable_name, cvar_cost_variable)
297      cvar_cost_variable.construct()
298
299      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
300      cvar_eta_variable = Var(name=cvar_eta_variable_name)
301      setattr(master_binding_instance, cvar_eta_variable_name, cvar_eta_variable)     
302      cvar_eta_variable.construct()
303
304   master_binding_instance.presolve()
305
306   # ditto above for the (non-expected) cost variable.
307   for stage in scenario_tree._stages:
308
309      (cost_variable,cost_variable_index) = stage._cost_variable
310
311      print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index     
312
313      for tree_node in stage._tree_nodes:
314
315         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
316
317         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)
318
319         # this is undoubtedly wasteful, in that a cost variable
320         # for each tree node is created with *all* indices.
321         new_cost_variable = copy.deepcopy(cost_variable)
322         new_cost_variable.name = new_cost_variable_name
323         new_cost_variable._model = master_binding_instance
324         setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)
325
326         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
327         new_cost_variable[cost_variable_index].var = new_cost_variable
328         if cost_variable_index is not None:
329            # if the variable index is None, the variable is derived from a VarValue, so the
330            # name gets updated automagically.
331            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
332
333         for scenario in tree_node._scenarios:
334
335            scenario_instance = instances[scenario._name]
336            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
337            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
338            new_constraint = Constraint(name=new_constraint_name)
339            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
340            new_constraint.add(None, (0.0, new_expr, 0.0))
341            new_constraint._model = master_binding_instance
342            setattr(master_binding_instance, new_constraint_name, new_constraint)
343
344   # create the constraints for computing the master per-node cost variables,
345   # i.e., the current node cost and the expected cost of the child nodes.
346   # if the root, then the constraint is just the objective.
347
348   for stage in scenario_tree._stages:
349
350      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
351
352      for tree_node in stage._tree_nodes:
353
354         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
355         node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
356
357         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
358         node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                       
359           
360         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
361
362         for child_node in tree_node._children:
363
364            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
365            child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
366            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
367
368         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
369         new_constraint = Constraint(name=new_constraint_name)
370         new_constraint.add(None, (0.0, constraint_expr, 0.0))
371         new_constraint._model = master_binding_instance                     
372         setattr(master_binding_instance, new_constraint_name, new_constraint)
373
374         if tree_node._parent is None:
375
376            an_instance = instances[instances.keys()[0]]
377            an_objective = an_instance.active_components(Objective)
378            opt_sense = an_objective[an_objective.keys()[0]].sense
379
380            opt_expression = node_expected_cost_variable
381
382            if generate_weighted_cvar is True:
383               cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
384               cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
385               opt_expression += cvar_weight * cvar_cost_variable
386
387            new_objective = Objective(name="MASTER", sense=opt_sense)
388            new_objective._data[None].expr = opt_expression
389            setattr(master_binding_instance, "MASTER", new_objective)
390
391   # CVaR requires the addition of a variable per scenario to represent the cost excess,
392   # and a constraint to compute the cost excess relative to eta. we also replicate (following
393   # what we do for node cost variables) an eta variable for each scenario instance, and
394   # require equality with the master eta variable via constraints.
395   if generate_weighted_cvar is True:
396     
397      root_node = scenario_tree._stages[0]._tree_nodes[0]
398
399      master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
400      master_cvar_eta_variable = getattr(master_binding_instance, master_cvar_eta_variable_name)
401     
402      for scenario_name in instances.keys():
403         scenario_instance = instances[scenario_name]
404
405         # unique names are required because the presolve isn't
406         # aware of the "owning" models for variables.
407         cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
408         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
409         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
410         cvar_excess_variable.construct()
411
412         cvar_eta_variable_name = "CVAR_ETA"
413         cvar_eta_variable = Var(name=cvar_eta_variable_name)
414         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
415         cvar_eta_variable.construct()
416
417         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
418         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
419         compute_excess_expression = cvar_excess_variable
420         for node in scenario_tree._scenario_map[scenario_name]._node_list:
421            (cost_variable, cost_variable_idx) = node._stage._cost_variable
422            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
423         compute_excess_expression += cvar_eta_variable
424         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
425         compute_excess_constraint._model = scenario_instance
426         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
427
428         eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
429         eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
430         eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
431         eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
432         eta_equality_constraint._model = master_binding_instance
433         setattr(master_binding_instance, eta_equality_constraint_name, eta_equality_constraint)
434
435      # add the constraint to compute the master CVaR variable value. iterate
436      # over scenario instances to create the expected excess component first.
437      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
438      cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
439      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
440      cvar_eta_variable = getattr(master_binding_instance, cvar_eta_variable_name)
441     
442      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
443     
444      for scenario_name in instances.keys():
445         scenario_instance = instances[scenario_name]
446         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
447
448         scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
449         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
450
451         cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
452
453      compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
454      compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
455      compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
456      compute_cvar_cost_constraint._model = master_binding_instance
457      setattr(master_binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
458
459   # after mucking with instances, presolve to collect terms required prior to output.         
460   # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios.
461   for scenario_name in instances.keys():
462      scenario_instance = instances[scenario_name]   
463      scenario_instance.presolve()         
464
465   master_binding_instance.presolve()
466
467   ################################################################################################
468   #### WRITE THE COMPOSITE MODEL #################################################################
469   ################################################################################################
470
471   print ""
472   print "Starting to write extensive form"
473
474   # create the output file.
475   problem_writer = cpxlp.ProblemWriter_cpxlp()
476   output_file = open(output_filename,"w")
477
478   problem_writer._output_prefixes = True # we always want prefixes
479
480   ################################################################################################
481   #### WRITE THE MASTER OBJECTIVE ################################################################
482   ################################################################################################
483
484   # write the objective for the master binding instance.
485   problem_writer._output_objectives = True
486   problem_writer._output_constraints = False
487   problem_writer._output_variables = False
488
489   print >>output_file, "\\ Begin objective block for master"
490   problem_writer._print_model_LP(master_binding_instance, output_file)
491   print >>output_file, "\\ End objective block for master"
492   print >>output_file, ""
493
494   ################################################################################################
495   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
496   ################################################################################################
497
498   print >>output_file, "s.t."
499   print >>output_file, ""
500   
501   problem_writer._output_objectives = False
502   problem_writer._output_constraints = True
503   problem_writer._output_variables = False
504
505   print >>output_file, "\\ Begin constraint block for master"
506   problem_writer._print_model_LP(master_binding_instance, output_file)
507   print >>output_file, "\\ End constraint block for master",
508   print >>output_file, ""
509
510   for scenario_name in instances.keys():
511      instance = instances[scenario_name]
512      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
513      problem_writer._print_model_LP(instance, output_file)
514      print >>output_file, "\\ End constraint block for scenario",scenario_name
515      print >>output_file, ""
516
517   ################################################################################################
518   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
519   ################################################################################################
520
521   # write the variables for the master binding instance, and then for each scenario.
522   print >>output_file, "bounds"
523   print >>output_file, ""
524   
525   problem_writer._output_objectives = False
526   problem_writer._output_constraints = False
527   problem_writer._output_variables = True
528
529   # first step: write variable bounds
530
531   problem_writer._output_continuous_variables = True
532   problem_writer._output_integer_variables = False
533   problem_writer._output_binary_variables = False
534
535   print >>output_file, "\\ Begin variable bounds block for master"
536   problem_writer._print_model_LP(master_binding_instance, output_file)
537   print >>output_file, "\\ End variable bounds block for master"
538   print >>output_file, ""
539
540   for scenario_name in instances.keys():
541      instance = instances[scenario_name]
542      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
543      problem_writer._print_model_LP(instance, output_file)
544      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
545      print >>output_file, ""
546
547   # second step: write integer indicators.
548
549   problem_writer._output_continuous_variables = False
550   problem_writer._output_integer_variables = True
551
552   if integers_present(master_binding_instance, instances) is True:
553
554      print >>output_file, "integer"
555      print >>output_file, ""
556
557      print >>output_file, "\\ Begin integer variable block for master"
558      problem_writer._print_model_LP(master_binding_instance, output_file)
559      print >>output_file, "\\ End integer variable block for master"
560      print >>output_file, ""
561   
562      for scenario_name in instances.keys():
563         instance = instances[scenario_name]
564         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name
565         problem_writer._print_model_LP(instance, output_file)
566         print >>output_file, "\\ End integer variable block for scenario",scenario_name
567         print >>output_file, ""
568
569   # third step: write binary indicators.
570
571   problem_writer._output_integer_variables = False
572   problem_writer._output_binary_variables = True
573
574   if binaries_present(master_binding_instance, instances) is True:
575
576      print >>output_file, "binary"
577      print >>output_file, ""
578
579      print >>output_file, "\\ Begin binary variable block for master"
580      problem_writer._print_model_LP(master_binding_instance, output_file)
581      print >>output_file, "\\ End binary variable block for master"
582      print >>output_file, ""
583   
584      for scenario_name in instances.keys():
585         instance = instances[scenario_name]
586         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
587         problem_writer._print_model_LP(instance, output_file)
588         print >>output_file, "\\ End integer binary block for scenario",scenario_name
589         print >>output_file, ""
590
591   # wrap up.
592   print >>output_file, "end"
593
594   # clean up.
595   output_file.close()
596
597   print ""
598   print "Output file written to file=",output_filename
599
600   print ""
601   print "Done..."
602
603   end_time = time.time()
604
605   print ""
606   print "Total execution time=%8.2f seconds" %(end_time - start_time)
607   print ""
608
609def write_ef(scenario_tree, instances, output_filename):
610
611   start_time = time.time()
612
613   ################################################################################################
614   #### CREATE THE MASTER / BINDING INSTANCE ######################################################
615   ################################################################################################
616
617   master_binding_instance = Model()
618   master_binding_instance.name = "MASTER"
619
620   # walk the scenario tree - create variables representing the common values for all scenarios
621   # associated with that node. the constraints will be created later. also create expected-cost
622   # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
623   # master variables are created for all nodes but those in the last stage. expected cost variables
624   # are, for no particularly good reason other than easy coding, created for nodes in all stages.
625   print "Creating variables for master binding instance"
626
627   for stage in scenario_tree._stages:
628
629      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
630
631         print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=", index_template
632
633         for tree_node in stage._tree_nodes:
634
635            if stage != scenario_tree._stages[-1]:     
636
637               master_variable_name = tree_node._name + "_" + stage_variable.name
638
639               # because there may be a single stage variable and multiple indices, check
640               # for the existence of the variable at this node - if you don't, you'll
641               # inadvertently over-write what was there previously!
642               master_variable = None
643               try:
644                  master_variable = getattr(master_binding_instance, master_variable_name)
645               except:
646                  # the deepcopy is probably too expensive (and unnecessary) computationally -
647                  # easier to just use the constructor with the stage variable index/bounds/etc.
648                  # NOTE: need to re-assign the master variables for each _varval - they probably
649                  #       point to a bogus model.
650                  new_master_variable = copy.deepcopy(stage_variable)
651                  new_master_variable.name = master_variable_name
652                  new_master_variable._model = master_binding_instance
653                  setattr(master_binding_instance, master_variable_name, new_master_variable)
654
655                  master_variable = new_master_variable
656
657               for index in stage_variable_indices:
658
659                  is_used = True # until proven otherwise                     
660                  for scenario in tree_node._scenarios:
661                     instance = instances[scenario._name]
662                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
663                        is_used = False
664
665                  is_fixed = False # until proven otherwise
666                  for scenario in tree_node._scenarios:
667                     instance = instances[scenario._name]
668                     if getattr(instance,stage_variable.name)[index].fixed is True:
669                        is_fixed = True
670
671                  if (is_used is True) and (is_fixed is False):
672                           
673                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
674                     # and because presolve/simplification is name-based, the names *have* to be different.
675                     master_variable[index].var = master_variable
676                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
677
678                     for scenario in tree_node._scenarios:
679
680                        scenario_instance = instances[scenario._name]
681                        scenario_variable = getattr(scenario_instance, stage_variable.name)
682                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
683                        new_constraint = Constraint(name=new_constraint_name)
684                        new_expr = master_variable[index] - scenario_variable[index]
685                        new_constraint.add(None, (0.0, new_expr, 0.0))
686                        new_constraint._model = master_binding_instance
687                        setattr(master_binding_instance, new_constraint_name, new_constraint)
688
689            # create a variable to represent the expected cost at this node -
690            # the constraint to compute this comes later.
691            expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
692            expected_cost_variable = Var(name=expected_cost_variable_name)
693            expected_cost_variable._model = master_binding_instance
694            setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
695
696   master_binding_instance.presolve()
697
698   # ditto above for the (non-expected) cost variable.
699   for stage in scenario_tree._stages:
700
701      (cost_variable,cost_variable_index) = stage._cost_variable
702
703      print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index     
704
705      for tree_node in stage._tree_nodes:
706
707         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
708
709         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)
710
711         # this is undoubtedly wasteful, in that a cost variable
712         # for each tree node is created with *all* indices.
713         new_cost_variable = copy.deepcopy(cost_variable)
714         new_cost_variable.name = new_cost_variable_name
715         new_cost_variable._model = master_binding_instance
716         setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)
717
718         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
719         new_cost_variable[cost_variable_index].var = new_cost_variable
720         if cost_variable_index is not None:
721            # if the variable index is None, the variable is derived from a VarValue, so the
722            # name gets updated automagically.
723            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
724
725         for scenario in tree_node._scenarios:
726
727            scenario_instance = instances[scenario._name]
728            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
729            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
730            new_constraint = Constraint(name=new_constraint_name)
731            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
732            new_constraint.add(None, (0.0, new_expr, 0.0))
733            new_constraint._model = master_binding_instance
734            setattr(master_binding_instance, new_constraint_name, new_constraint)
735
736   # create the constraints for computing the master per-node cost variables,
737   # i.e., the current node cost and the expected cost of the child nodes.
738   # if the root, then the constraint is just the objective.
739
740   for stage in scenario_tree._stages:
741
742      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
743
744      for tree_node in stage._tree_nodes:
745
746         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
747         node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
748
749         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
750         node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                       
751           
752         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
753
754         for child_node in tree_node._children:
755
756            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
757            child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
758            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
759
760         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
761         new_constraint = Constraint(name=new_constraint_name)
762         new_constraint.add(None, (0.0, constraint_expr, 0.0))
763         new_constraint._model = master_binding_instance                     
764         setattr(master_binding_instance, new_constraint_name, new_constraint)
765
766         if tree_node._parent is None:
767
768            an_instance = instances[instances.keys()[0]]
769            an_objective = an_instance.active_components(Objective)
770            opt_sense = an_objective[an_objective.keys()[0]].sense
771
772            new_objective = Objective(name="MASTER", sense=opt_sense)
773            new_objective._data[None].expr = node_expected_cost_variable
774            setattr(master_binding_instance, "MASTER", new_objective)
775
776   master_binding_instance.presolve()
777
778   ################################################################################################
779   #### WRITE THE COMPOSITE MODEL #################################################################
780   ################################################################################################
781
782   print ""
783   print "Starting to write extensive form"
784
785   # create the output file.
786   problem_writer = cpxlp.ProblemWriter_cpxlp()
787   output_file = open(output_filename,"w")
788
789   problem_writer._output_prefixes = True # we always want prefixes
790
791   ################################################################################################
792   #### WRITE THE MASTER OBJECTIVE ################################################################
793   ################################################################################################
794
795   # write the objective for the master binding instance.
796   problem_writer._output_objectives = True
797   problem_writer._output_constraints = False
798   problem_writer._output_variables = False
799
800   print >>output_file, "\\ Begin objective block for master"
801   problem_writer._print_model_LP(master_binding_instance, output_file)
802   print >>output_file, "\\ End objective block for master"
803   print >>output_file, ""
804
805   ################################################################################################
806   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
807   ################################################################################################
808
809   print >>output_file, "s.t."
810   print >>output_file, ""
811   
812   problem_writer._output_objectives = False
813   problem_writer._output_constraints = True
814   problem_writer._output_variables = False
815
816   print >>output_file, "\\ Begin constraint block for master"
817   problem_writer._print_model_LP(master_binding_instance, output_file)
818   print >>output_file, "\\ End constraint block for master",
819   print >>output_file, ""
820
821   for scenario_name in instances.keys():
822      instance = instances[scenario_name]
823      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
824      problem_writer._print_model_LP(instance, output_file)
825      print >>output_file, "\\ End constraint block for scenario",scenario_name
826      print >>output_file, ""
827
828   ################################################################################################
829   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
830   ################################################################################################
831
832   # write the variables for the master binding instance, and then for each scenario.
833   print >>output_file, "bounds"
834   print >>output_file, ""
835   
836   problem_writer._output_objectives = False
837   problem_writer._output_constraints = False
838   problem_writer._output_variables = True
839
840   # first step: write variable bounds
841
842   problem_writer._output_continuous_variables = True
843   problem_writer._output_integer_variables = False
844   problem_writer._output_binary_variables = False
845
846   print >>output_file, "\\ Begin variable bounds block for master"
847   problem_writer._print_model_LP(master_binding_instance, output_file)
848   print >>output_file, "\\ End variable bounds block for master"
849   print >>output_file, ""
850   
851   for scenario_name in instances.keys():
852      instance = instances[scenario_name]
853      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
854      problem_writer._print_model_LP(instance, output_file)
855      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
856      print >>output_file, ""
857
858   # second step: write integer indicators.
859
860   problem_writer._output_continuous_variables = False
861   problem_writer._output_integer_variables = True
862
863   if integers_present(master_binding_instance, instances) is True:
864
865      print >>output_file, "integer"
866      print >>output_file, ""
867
868      print >>output_file, "\\ Begin integer variable block for master"
869      problem_writer._print_model_LP(master_binding_instance, output_file)
870      print >>output_file, "\\ End integer variable block for master"
871      print >>output_file, ""
872   
873      for scenario_name in instances.keys():
874         instance = instances[scenario_name]
875         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name
876         problem_writer._print_model_LP(instance, output_file)
877         print >>output_file, "\\ End integer variable block for scenario",scenario_name
878         print >>output_file, ""
879
880   # third step: write binary indicators.
881
882   problem_writer._output_integer_variables = False
883   problem_writer._output_binary_variables = True
884
885   if binaries_present(master_binding_instance, instances) is True:
886
887      print >>output_file, "binary"
888      print >>output_file, ""
889
890      print >>output_file, "\\ Begin binary variable block for master"
891      problem_writer._print_model_LP(master_binding_instance, output_file)
892      print >>output_file, "\\ End binary variable block for master"
893      print >>output_file, ""
894   
895      for scenario_name in instances.keys():
896         instance = instances[scenario_name]
897         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
898         problem_writer._print_model_LP(instance, output_file)
899         print >>output_file, "\\ End integer binary block for scenario",scenario_name
900         print >>output_file, ""
901
902   # wrap up.
903   print >>output_file, "end"
904
905   # clean up.
906   output_file.close()
907
908   print ""
909   print "Output file written to file=",output_filename
910
911   print ""
912   print "Done..."
913
914   end_time = time.time()
915
916   print ""
917   print "Total execution time=%8.2f seconds" %(end_time - start_time)
918   print ""   
Note: See TracBrowser for help on using the repository browser.