source: coopr.pysp/stable/2.3/coopr/pysp/ef.py @ 2317

Last change on this file since 2317 was 2317, checked in by wehart, 10 years ago

Merged revisions 2246-2316 via svnmerge from
https://software.sandia.gov/svn/public/coopr/coopr.pysp/trunk

........

r2246 | wehart | 2010-02-01 21:10:48 -0700 (Mon, 01 Feb 2010) | 2 lines


Tagging coopr.pysp release 2.2

........

r2247 | wehart | 2010-02-01 21:46:08 -0700 (Mon, 01 Feb 2010) | 2 lines


Documentation fix.

........

r2248 | jwatson | 2010-02-03 08:44:59 -0700 (Wed, 03 Feb 2010) | 3 lines


Adding os.path.expanduser wrappers around all directory/filenames, to facilitate correct processing of ~ characters.

........

r2249 | jwatson | 2010-02-03 09:28:25 -0700 (Wed, 03 Feb 2010) | 3 lines


Misc fixes.

........

r2254 | jwatson | 2010-02-03 21:09:38 -0700 (Wed, 03 Feb 2010) | 3 lines


Changed PyomoModelData? call from add_data_mumble() to add().

........

r2255 | jwatson | 2010-02-03 21:38:35 -0700 (Wed, 03 Feb 2010) | 3 lines


Re-factoring of PH options parser code to accomodate MRP work being done by DLW.

........

r2256 | jwatson | 2010-02-03 22:18:11 -0700 (Wed, 03 Feb 2010) | 3 lines


Major speed improvements in the EF writer by avoiding Python deep-copes - saves a few orders of magnitude of run-time.

........

r2260 | jwatson | 2010-02-04 21:44:29 -0700 (Thu, 04 Feb 2010) | 3 lines


Refactoring of ph initialization routines to support sampling and bundling.

........

r2261 | jwatson | 2010-02-04 21:46:37 -0700 (Thu, 04 Feb 2010) | 3 lines


Renaming ph_script module to phinit, which is more accurate with the newly factored, library-like functionality.

........

r2262 | jwatson | 2010-02-04 21:54:50 -0700 (Thu, 04 Feb 2010) | 3 lines


Minor improvement to phinit functionality.

........

r2267 | jwatson | 2010-02-05 15:09:02 -0700 (Fri, 05 Feb 2010) | 3 lines


Initial PySP unit tests!!!

........

r2268 | wehart | 2010-02-05 15:25:52 -0700 (Fri, 05 Feb 2010) | 2 lines


A fix to the tests.

........

r2269 | jwatson | 2010-02-05 15:46:35 -0700 (Fri, 05 Feb 2010) | 3 lines


Added SIZES3 PySP test and added absolute paths to "runph" script.

........

r2270 | jwatson | 2010-02-05 16:01:08 -0700 (Fri, 05 Feb 2010) | 3 lines


PySP tests need absolute output paths!

........

r2271 | jwatson | 2010-02-05 16:20:01 -0700 (Fri, 05 Feb 2010) | 3 lines


Making unit tests for PySP compatible with coverage utilities. Farmer examples work, SIZES3 doesn't for some reason.

........

r2272 | jwatson | 2010-02-05 16:59:46 -0700 (Fri, 05 Feb 2010) | 3 lines


Various fixes to PySP unit tests. Changing name of testphextension to examplephextension - with the test prefix, coverage tests import the module, which causes all kinds of issues.

........

r2275 | jwatson | 2010-02-06 13:45:11 -0700 (Sat, 06 Feb 2010) | 3 lines


Testing improvements. From lpython, the tests run individually just fine. In aggregate, only the first run passes - independent of what test that might be! Something to stare at later....

........

r2277 | jwatson | 2010-02-06 23:30:13 -0700 (Sat, 06 Feb 2010) | 3 lines


Fixed problem with runph --profile option, broken by my recent factoring of phinit.py.

........

r2288 | jwatson | 2010-02-08 13:32:41 -0700 (Mon, 08 Feb 2010) | 3 lines


Significant initialization speed reductions in the WW PH extension for PySP.

........

r2290 | jwatson | 2010-02-08 19:22:32 -0700 (Mon, 08 Feb 2010) | 1 line


Initial commit of multi-stage capacity expansion problem in PySP

........

r2291 | jwatson | 2010-02-08 19:23:24 -0700 (Mon, 08 Feb 2010) | 1 line


Miscellaneous fix to ef writer involving indexed cost variables

........

r2296 | jwatson | 2010-02-09 18:56:48 -0700 (Tue, 09 Feb 2010) | 3 lines


Removing monster-sized LP file from the PySP forestry examples directory.

........

r2297 | jwatson | 2010-02-09 18:59:05 -0700 (Tue, 09 Feb 2010) | 3 lines


The extensive forms in the PySP forestry example were massive - and are now gone.

........

r2298 | jwatson | 2010-02-09 19:03:09 -0700 (Tue, 09 Feb 2010) | 1 line


Removing output logs for PySP network flow example

........

r2299 | jwatson | 2010-02-09 19:04:12 -0700 (Tue, 09 Feb 2010) | 1 line


Removing a big network flow EF

........

r2300 | jwatson | 2010-02-09 19:05:08 -0700 (Tue, 09 Feb 2010) | 3 lines


Removing PySP cap example EF to free up space.

........

r2301 | jwatson | 2010-02-09 19:14:26 -0700 (Tue, 09 Feb 2010) | 1 line


Performance improvements to PH obtained by processing scenario sub-problem results as they come in, instead of waiting for them after a solver barrier sync

........

File size: 42.6 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(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                  new_master_variable_index = stage_variable._index
240                  new_master_variable = None
241                  if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
242                     new_master_variable = Var(name=stage_variable.name)
243                  else:
244                     new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
245                  new_master_variable.construct()
246                  new_master_variable._model = master_binding_instance
247                  setattr(master_binding_instance, master_variable_name, new_master_variable)
248
249                  # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS!
250
251                  master_variable = new_master_variable
252
253               for index in stage_variable_indices:
254
255                  is_used = True # until proven otherwise                     
256                  for scenario in tree_node._scenarios:
257                     instance = instances[scenario._name]
258                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
259                        is_used = False
260
261                  is_fixed = False # until proven otherwise
262                  for scenario in tree_node._scenarios:
263                     instance = instances[scenario._name]
264                     if getattr(instance,stage_variable.name)[index].fixed is True:
265                        is_fixed = True
266
267                  if (is_used is True) and (is_fixed is False):
268                           
269                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
270                     # and because presolve/simplification is name-based, the names *have* to be different.
271                     master_variable[index].var = master_variable
272                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
273
274                     for scenario in tree_node._scenarios:
275
276                        scenario_instance = instances[scenario._name]
277                        scenario_variable = getattr(scenario_instance, stage_variable.name)
278                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
279                        new_constraint = Constraint(name=new_constraint_name)
280                        new_expr = master_variable[index] - scenario_variable[index]
281                        new_constraint.add(None, (0.0, new_expr, 0.0))
282                        new_constraint._model = master_binding_instance
283                        setattr(master_binding_instance, new_constraint_name, new_constraint)
284
285            # create a variable to represent the expected cost at this node -
286            # the constraint to compute this comes later.
287            expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
288            expected_cost_variable = Var(name=expected_cost_variable_name)
289            expected_cost_variable._model = master_binding_instance
290            setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
291
292   # if we're generating the weighted CVaR objective term, create the corresponding variable and
293   # the master CVaR eta variable.
294   if generate_weighted_cvar is True:
295      root_node = scenario_tree._stages[0]._tree_nodes[0]
296     
297      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
298      cvar_cost_variable = Var(name=cvar_cost_variable_name)
299      setattr(master_binding_instance, cvar_cost_variable_name, cvar_cost_variable)
300      cvar_cost_variable.construct()
301
302      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
303      cvar_eta_variable = Var(name=cvar_eta_variable_name)
304      setattr(master_binding_instance, cvar_eta_variable_name, cvar_eta_variable)     
305      cvar_eta_variable.construct()
306
307   master_binding_instance.presolve()
308
309   # ditto above for the (non-expected) cost variable.
310   for stage in scenario_tree._stages:
311
312      (cost_variable,cost_variable_index) = stage._cost_variable
313
314      print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index     
315
316      for tree_node in stage._tree_nodes:
317
318         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)         
319
320         # this is undoubtedly wasteful, in that a cost variable
321         # for each tree node is created with *all* indices.         
322         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
323         new_cost_variable_index = cost_variable._index
324         new_cost_variable = None
325         if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
326            new_cost_variable = Var(name=new_cost_variable_name)
327         else:
328            new_cost_variable = Var(new_cost_variable_index, name=new_cost_variable_name)
329         new_cost_variable.construct()
330         new_cost_variable._model = master_binding_instance
331         setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)         
332
333         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
334         new_cost_variable[cost_variable_index].var = new_cost_variable
335         if cost_variable_index is not None:
336            # if the variable index is None, the variable is derived from a VarValue, so the
337            # name gets updated automagically.
338            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
339
340         for scenario in tree_node._scenarios:
341
342            scenario_instance = instances[scenario._name]
343            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
344            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
345            new_constraint = Constraint(name=new_constraint_name)
346            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
347            new_constraint.add(None, (0.0, new_expr, 0.0))
348            new_constraint._model = master_binding_instance
349            setattr(master_binding_instance, new_constraint_name, new_constraint)
350
351   # create the constraints for computing the master per-node cost variables,
352   # i.e., the current node cost and the expected cost of the child nodes.
353   # if the root, then the constraint is just the objective.
354
355   for stage in scenario_tree._stages:
356
357      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
358
359      for tree_node in stage._tree_nodes:
360
361         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
362         node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
363
364         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
365         node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                       
366           
367         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
368
369         for child_node in tree_node._children:
370
371            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
372            child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
373            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
374
375         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
376         new_constraint = Constraint(name=new_constraint_name)
377         new_constraint.add(None, (0.0, constraint_expr, 0.0))
378         new_constraint._model = master_binding_instance                     
379         setattr(master_binding_instance, new_constraint_name, new_constraint)
380
381         if tree_node._parent is None:
382
383            an_instance = instances[instances.keys()[0]]
384            an_objective = an_instance.active_components(Objective)
385            opt_sense = an_objective[an_objective.keys()[0]].sense
386
387            opt_expression = node_expected_cost_variable
388
389            if generate_weighted_cvar is True:
390               cvar_cost_variable_name = "CVAR_COST_" + tree_node._name
391               cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
392               opt_expression += cvar_weight * cvar_cost_variable
393
394            new_objective = Objective(name="MASTER", sense=opt_sense)
395            new_objective._data[None].expr = opt_expression
396            setattr(master_binding_instance, "MASTER", new_objective)
397
398   # CVaR requires the addition of a variable per scenario to represent the cost excess,
399   # and a constraint to compute the cost excess relative to eta. we also replicate (following
400   # what we do for node cost variables) an eta variable for each scenario instance, and
401   # require equality with the master eta variable via constraints.
402   if generate_weighted_cvar is True:
403     
404      root_node = scenario_tree._stages[0]._tree_nodes[0]
405
406      master_cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
407      master_cvar_eta_variable = getattr(master_binding_instance, master_cvar_eta_variable_name)
408     
409      for scenario_name in instances.keys():
410         scenario_instance = instances[scenario_name]
411
412         # unique names are required because the presolve isn't
413         # aware of the "owning" models for variables.
414         cvar_excess_variable_name = "CVAR_EXCESS_"+scenario_name
415         cvar_excess_variable = Var(name=cvar_excess_variable_name, domain=NonNegativeReals)
416         setattr(scenario_instance, cvar_excess_variable_name, cvar_excess_variable)
417         cvar_excess_variable.construct()
418
419         cvar_eta_variable_name = "CVAR_ETA"
420         cvar_eta_variable = Var(name=cvar_eta_variable_name)
421         setattr(scenario_instance, cvar_eta_variable_name, cvar_eta_variable)
422         cvar_eta_variable.construct()
423
424         compute_excess_constraint_name = "COMPUTE_SCENARIO_EXCESS"
425         compute_excess_constraint = Constraint(name=compute_excess_constraint_name)
426         compute_excess_expression = cvar_excess_variable
427         for node in scenario_tree._scenario_map[scenario_name]._node_list:
428            (cost_variable, cost_variable_idx) = node._stage._cost_variable
429            compute_excess_expression -= getattr(scenario_instance, cost_variable.name)[cost_variable_idx]
430         compute_excess_expression += cvar_eta_variable
431         compute_excess_constraint.add(None, (0.0, compute_excess_expression, None))
432         compute_excess_constraint._model = scenario_instance
433         setattr(scenario_instance, compute_excess_constraint_name, compute_excess_constraint)
434
435         eta_equality_constraint_name = "MASTER_ETA_EQUALITY_WITH_" + scenario_instance.name
436         eta_equality_constraint = Constraint(name=eta_equality_constraint_name)
437         eta_equality_expr = master_cvar_eta_variable - cvar_eta_variable
438         eta_equality_constraint.add(None, (0.0, eta_equality_expr, 0.0))
439         eta_equality_constraint._model = master_binding_instance
440         setattr(master_binding_instance, eta_equality_constraint_name, eta_equality_constraint)
441
442      # add the constraint to compute the master CVaR variable value. iterate
443      # over scenario instances to create the expected excess component first.
444      cvar_cost_variable_name = "CVAR_COST_" + root_node._name
445      cvar_cost_variable = getattr(master_binding_instance, cvar_cost_variable_name)
446      cvar_eta_variable_name = "CVAR_ETA_" + root_node._name
447      cvar_eta_variable = getattr(master_binding_instance, cvar_eta_variable_name)
448     
449      cvar_cost_expression = cvar_cost_variable - cvar_eta_variable
450     
451      for scenario_name in instances.keys():
452         scenario_instance = instances[scenario_name]
453         scenario_probability = scenario_tree._scenario_map[scenario_name]._probability
454
455         scenario_excess_variable_name = "CVAR_EXCESS_"+scenario_name
456         scenario_excess_variable = getattr(scenario_instance, scenario_excess_variable_name)
457
458         cvar_cost_expression = cvar_cost_expression - (scenario_probability * scenario_excess_variable) / (1.0 - risk_alpha)
459
460      compute_cvar_cost_constraint_name = "COMPUTE_CVAR_COST"
461      compute_cvar_cost_constraint = Constraint(name=compute_cvar_cost_constraint_name)
462      compute_cvar_cost_constraint.add(None, (0.0, cvar_cost_expression, 0.0))
463      compute_cvar_cost_constraint._model = master_binding_instance
464      setattr(master_binding_instance, compute_cvar_cost_constraint_name, compute_cvar_cost_constraint)
465
466   # after mucking with instances, presolve to collect terms required prior to output.         
467   # IMPT: Do the scenario instances first, as the master depends on variables in the scenarios.
468   for scenario_name in instances.keys():
469      scenario_instance = instances[scenario_name]   
470      scenario_instance.presolve()         
471
472   master_binding_instance.presolve()
473
474   ################################################################################################
475   #### WRITE THE COMPOSITE MODEL #################################################################
476   ################################################################################################
477
478   print ""
479   print "Starting to write extensive form"
480
481   # create the output file.
482   problem_writer = cpxlp.ProblemWriter_cpxlp()
483   output_file = open(output_filename,"w")
484
485   problem_writer._output_prefixes = True # we always want prefixes
486
487   ################################################################################################
488   #### WRITE THE MASTER OBJECTIVE ################################################################
489   ################################################################################################
490
491   # write the objective for the master binding instance.
492   problem_writer._output_objectives = True
493   problem_writer._output_constraints = False
494   problem_writer._output_variables = False
495
496   print >>output_file, "\\ Begin objective block for master"
497   problem_writer._print_model_LP(master_binding_instance, output_file)
498   print >>output_file, "\\ End objective block for master"
499   print >>output_file, ""
500
501   ################################################################################################
502   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
503   ################################################################################################
504
505   print >>output_file, "s.t."
506   print >>output_file, ""
507   
508   problem_writer._output_objectives = False
509   problem_writer._output_constraints = True
510   problem_writer._output_variables = False
511
512   print >>output_file, "\\ Begin constraint block for master"
513   problem_writer._print_model_LP(master_binding_instance, output_file)
514   print >>output_file, "\\ End constraint block for master",
515   print >>output_file, ""
516
517   for scenario_name in instances.keys():
518      instance = instances[scenario_name]
519      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
520      problem_writer._print_model_LP(instance, output_file)
521      print >>output_file, "\\ End constraint block for scenario",scenario_name
522      print >>output_file, ""
523
524   ################################################################################################
525   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
526   ################################################################################################
527
528   # write the variables for the master binding instance, and then for each scenario.
529   print >>output_file, "bounds"
530   print >>output_file, ""
531   
532   problem_writer._output_objectives = False
533   problem_writer._output_constraints = False
534   problem_writer._output_variables = True
535
536   # first step: write variable bounds
537
538   problem_writer._output_continuous_variables = True
539   problem_writer._output_integer_variables = False
540   problem_writer._output_binary_variables = False
541
542   print >>output_file, "\\ Begin variable bounds block for master"
543   problem_writer._print_model_LP(master_binding_instance, output_file)
544   print >>output_file, "\\ End variable bounds block for master"
545   print >>output_file, ""
546
547   for scenario_name in instances.keys():
548      instance = instances[scenario_name]
549      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
550      problem_writer._print_model_LP(instance, output_file)
551      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
552      print >>output_file, ""
553
554   # second step: write integer indicators.
555
556   problem_writer._output_continuous_variables = False
557   problem_writer._output_integer_variables = True
558
559   if integers_present(master_binding_instance, instances) is True:
560
561      print >>output_file, "integer"
562      print >>output_file, ""
563
564      print >>output_file, "\\ Begin integer variable block for master"
565      problem_writer._print_model_LP(master_binding_instance, output_file)
566      print >>output_file, "\\ End integer variable block for master"
567      print >>output_file, ""
568   
569      for scenario_name in instances.keys():
570         instance = instances[scenario_name]
571         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name
572         problem_writer._print_model_LP(instance, output_file)
573         print >>output_file, "\\ End integer variable block for scenario",scenario_name
574         print >>output_file, ""
575
576   # third step: write binary indicators.
577
578   problem_writer._output_integer_variables = False
579   problem_writer._output_binary_variables = True
580
581   if binaries_present(master_binding_instance, instances) is True:
582
583      print >>output_file, "binary"
584      print >>output_file, ""
585
586      print >>output_file, "\\ Begin binary variable block for master"
587      problem_writer._print_model_LP(master_binding_instance, output_file)
588      print >>output_file, "\\ End binary variable block for master"
589      print >>output_file, ""
590   
591      for scenario_name in instances.keys():
592         instance = instances[scenario_name]
593         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
594         problem_writer._print_model_LP(instance, output_file)
595         print >>output_file, "\\ End integer binary block for scenario",scenario_name
596         print >>output_file, ""
597
598   # wrap up.
599   print >>output_file, "end"
600
601   # clean up.
602   output_file.close()
603
604   print ""
605   print "Output file written to file=",output_filename
606
607   print ""
608   print "Done..."
609
610   end_time = time.time()
611
612   print ""
613   print "Total execution time=%8.2f seconds" %(end_time - start_time)
614   print ""
615
616def write_ef(scenario_tree, instances, output_filename):
617
618   start_time = time.time()
619
620   ################################################################################################
621   #### CREATE THE MASTER / BINDING INSTANCE ######################################################
622   ################################################################################################
623
624   master_binding_instance = Model()
625   master_binding_instance.name = "MASTER"
626
627   # walk the scenario tree - create variables representing the common values for all scenarios
628   # associated with that node. the constraints will be created later. also create expected-cost
629   # variables for each node, to be computed via constraints/objectives defined in a subsequent pass.
630   # master variables are created for all nodes but those in the last stage. expected cost variables
631   # are, for no particularly good reason other than easy coding, created for nodes in all stages.
632   print "Creating variables for master binding instance"
633
634   for stage in scenario_tree._stages:
635
636      for (stage_variable, index_template, stage_variable_indices) in stage._variables:
637
638         print "Creating master variable and blending constraints for decision variable=", stage_variable, ", indices=", index_template
639
640         for tree_node in stage._tree_nodes:
641
642            if stage != scenario_tree._stages[-1]:     
643
644               master_variable_name = tree_node._name + "_" + stage_variable.name
645
646               # because there may be a single stage variable and multiple indices, check
647               # for the existence of the variable at this node - if you don't, you'll
648               # inadvertently over-write what was there previously!
649               master_variable = None
650               try:
651                  master_variable = getattr(master_binding_instance, master_variable_name)
652               except:
653                  new_master_variable_index = stage_variable._index
654                  new_master_variable = None
655                  if (len(new_master_variable_index) is 1) and (None in new_master_variable_index):
656                     new_master_variable = Var(name=stage_variable.name)
657                  else:
658                     new_master_variable = Var(new_master_variable_index, name=stage_variable.name)
659                  new_master_variable.construct()
660                  new_master_variable._model = master_binding_instance
661                  setattr(master_binding_instance, master_variable_name, new_master_variable)
662
663                  # TBD - TECHNICALLY, WE NEED TO COPY BOUNDS - BUT WE REALLY DON'T, AS THEY ARE ON THE PER-INSTNACE VARS!
664
665                  master_variable = new_master_variable
666
667               for index in stage_variable_indices:
668
669                  is_used = True # until proven otherwise                     
670                  for scenario in tree_node._scenarios:
671                     instance = instances[scenario._name]
672                     if getattr(instance,stage_variable.name)[index].status == VarStatus.unused:
673                        is_used = False
674
675                  is_fixed = False # until proven otherwise
676                  for scenario in tree_node._scenarios:
677                     instance = instances[scenario._name]
678                     if getattr(instance,stage_variable.name)[index].fixed is True:
679                        is_fixed = True
680
681                  if (is_used is True) and (is_fixed is False):
682                           
683                     # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
684                     # and because presolve/simplification is name-based, the names *have* to be different.
685                     master_variable[index].var = master_variable
686                     master_variable[index].name = tree_node._name + "_" + master_variable[index].name
687
688                     for scenario in tree_node._scenarios:
689
690                        scenario_instance = instances[scenario._name]
691                        scenario_variable = getattr(scenario_instance, stage_variable.name)
692                        new_constraint_name = scenario._name + "_" + master_variable_name + "_" + str(index)
693                        new_constraint = Constraint(name=new_constraint_name)
694                        new_expr = master_variable[index] - scenario_variable[index]
695                        new_constraint.add(None, (0.0, new_expr, 0.0))
696                        new_constraint._model = master_binding_instance
697                        setattr(master_binding_instance, new_constraint_name, new_constraint)
698
699            # create a variable to represent the expected cost at this node -
700            # the constraint to compute this comes later.
701            expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
702            expected_cost_variable = Var(name=expected_cost_variable_name)
703            expected_cost_variable._model = master_binding_instance
704            setattr(master_binding_instance, expected_cost_variable_name, expected_cost_variable)
705
706   master_binding_instance.presolve()
707
708   # ditto above for the (non-expected) cost variable.
709   for stage in scenario_tree._stages:
710
711      (cost_variable,cost_variable_index) = stage._cost_variable
712
713      print "Creating master variable and blending constraints for cost variable=", cost_variable, ", index=", cost_variable_index     
714
715      for tree_node in stage._tree_nodes:
716
717         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
718
719         # TBD - the following is bad - check to see if it's already there (I suspect some of them are!!!)
720
721         # this is undoubtedly wasteful, in that a cost variable
722         # for each tree node is created with *all* indices.
723         new_cost_variable_name = tree_node._name + "_" + cost_variable.name
724         new_cost_variable_index = cost_variable._index
725         new_cost_variable = None
726         if (len(new_cost_variable_index) is 1) and (None in new_cost_variable_index):
727            new_cost_variable = Var(name=new_cost_variable_name)
728         else:
729            new_cost_variable = Var(new_cost_variable_index, new_cost_variable_name)
730         new_cost_variable.construct()
731         new_cost_variable._model = master_binding_instance
732         setattr(master_binding_instance, new_cost_variable_name, new_cost_variable)                 
733
734         # the following is necessary, specifically to get the name - deepcopy won't reset these attributes.
735         new_cost_variable[cost_variable_index].var = new_cost_variable
736         if cost_variable_index is not None:
737            # if the variable index is None, the variable is derived from a VarValue, so the
738            # name gets updated automagically.
739            new_cost_variable[cost_variable_index].name = tree_node._name + "_" + new_cost_variable[cost_variable_index].name
740
741         for scenario in tree_node._scenarios:
742
743            scenario_instance = instances[scenario._name]
744            scenario_cost_variable = getattr(scenario_instance, cost_variable.name)
745            new_constraint_name = scenario._name + "_" + new_cost_variable_name + "_" + str(cost_variable_index)
746            new_constraint = Constraint(name=new_constraint_name)
747            new_expr = new_cost_variable[cost_variable_index] - scenario_cost_variable[cost_variable_index]
748            new_constraint.add(None, (0.0, new_expr, 0.0))
749            new_constraint._model = master_binding_instance
750            setattr(master_binding_instance, new_constraint_name, new_constraint)
751
752   # create the constraints for computing the master per-node cost variables,
753   # i.e., the current node cost and the expected cost of the child nodes.
754   # if the root, then the constraint is just the objective.
755
756   for stage in scenario_tree._stages:
757
758      (stage_cost_variable,stage_cost_variable_index) = stage._cost_variable
759
760      for tree_node in stage._tree_nodes:
761
762         node_expected_cost_variable_name = "EXPECTED_COST_" + tree_node._name
763         node_expected_cost_variable = getattr(master_binding_instance, node_expected_cost_variable_name)
764
765         node_cost_variable_name = tree_node._name + "_" + stage_cost_variable.name
766         node_cost_variable = getattr(master_binding_instance, node_cost_variable_name)                       
767           
768         constraint_expr = node_expected_cost_variable - node_cost_variable[stage_cost_variable_index]
769
770         for child_node in tree_node._children:
771
772            child_node_expected_cost_variable_name = "EXPECTED_COST_" + child_node._name
773            child_node_expected_cost_variable = getattr(master_binding_instance, child_node_expected_cost_variable_name)
774            constraint_expr = constraint_expr - (child_node._conditional_probability * child_node_expected_cost_variable)
775
776         new_constraint_name = "COST" + "_" + node_cost_variable_name + "_" + str(cost_variable_index)
777         new_constraint = Constraint(name=new_constraint_name)
778         new_constraint.add(None, (0.0, constraint_expr, 0.0))
779         new_constraint._model = master_binding_instance                     
780         setattr(master_binding_instance, new_constraint_name, new_constraint)
781
782         if tree_node._parent is None:
783
784            an_instance = instances[instances.keys()[0]]
785            an_objective = an_instance.active_components(Objective)
786            opt_sense = an_objective[an_objective.keys()[0]].sense
787
788            new_objective = Objective(name="MASTER", sense=opt_sense)
789            new_objective._data[None].expr = node_expected_cost_variable
790            setattr(master_binding_instance, "MASTER", new_objective)
791
792   master_binding_instance.presolve()
793
794   ################################################################################################
795   #### WRITE THE COMPOSITE MODEL #################################################################
796   ################################################################################################
797
798   print ""
799   print "Starting to write extensive form"
800
801   # create the output file.
802   problem_writer = cpxlp.ProblemWriter_cpxlp()
803   output_file = open(output_filename,"w")
804
805   problem_writer._output_prefixes = True # we always want prefixes
806
807   ################################################################################################
808   #### WRITE THE MASTER OBJECTIVE ################################################################
809   ################################################################################################
810
811   # write the objective for the master binding instance.
812   problem_writer._output_objectives = True
813   problem_writer._output_constraints = False
814   problem_writer._output_variables = False
815
816   print >>output_file, "\\ Begin objective block for master"
817   problem_writer._print_model_LP(master_binding_instance, output_file)
818   print >>output_file, "\\ End objective block for master"
819   print >>output_file, ""
820
821   ################################################################################################
822   #### WRITE THE CONSTRAINTS FOR THE MASTER MODEL AND ALL SCENARIO MODELS ########################
823   ################################################################################################
824
825   print >>output_file, "s.t."
826   print >>output_file, ""
827   
828   problem_writer._output_objectives = False
829   problem_writer._output_constraints = True
830   problem_writer._output_variables = False
831
832   print >>output_file, "\\ Begin constraint block for master"
833   problem_writer._print_model_LP(master_binding_instance, output_file)
834   print >>output_file, "\\ End constraint block for master",
835   print >>output_file, ""
836
837   for scenario_name in instances.keys():
838      instance = instances[scenario_name]
839      print >>output_file, "\\ Begin constraint block for scenario",scenario_name       
840      problem_writer._print_model_LP(instance, output_file)
841      print >>output_file, "\\ End constraint block for scenario",scenario_name
842      print >>output_file, ""
843
844   ################################################################################################
845   #### WRITE THE VARIABLES FOR THE MASTER MODEL AND ALL SCENARIO MODELS ##########################
846   ################################################################################################
847
848   # write the variables for the master binding instance, and then for each scenario.
849   print >>output_file, "bounds"
850   print >>output_file, ""
851   
852   problem_writer._output_objectives = False
853   problem_writer._output_constraints = False
854   problem_writer._output_variables = True
855
856   # first step: write variable bounds
857
858   problem_writer._output_continuous_variables = True
859   problem_writer._output_integer_variables = False
860   problem_writer._output_binary_variables = False
861
862   print >>output_file, "\\ Begin variable bounds block for master"
863   problem_writer._print_model_LP(master_binding_instance, output_file)
864   print >>output_file, "\\ End variable bounds block for master"
865   print >>output_file, ""
866   
867   for scenario_name in instances.keys():
868      instance = instances[scenario_name]
869      print >>output_file, "\\ Begin variable bounds block for scenario",scenario_name
870      problem_writer._print_model_LP(instance, output_file)
871      print >>output_file, "\\ End variable bounds block for scenario",scenario_name
872      print >>output_file, ""
873
874   # second step: write integer indicators.
875
876   problem_writer._output_continuous_variables = False
877   problem_writer._output_integer_variables = True
878
879   if integers_present(master_binding_instance, instances) is True:
880
881      print >>output_file, "integer"
882      print >>output_file, ""
883
884      print >>output_file, "\\ Begin integer variable block for master"
885      problem_writer._print_model_LP(master_binding_instance, output_file)
886      print >>output_file, "\\ End integer variable block for master"
887      print >>output_file, ""
888   
889      for scenario_name in instances.keys():
890         instance = instances[scenario_name]
891         print >>output_file, "\\ Begin integer variable block for scenario",scenario_name
892         problem_writer._print_model_LP(instance, output_file)
893         print >>output_file, "\\ End integer variable block for scenario",scenario_name
894         print >>output_file, ""
895
896   # third step: write binary indicators.
897
898   problem_writer._output_integer_variables = False
899   problem_writer._output_binary_variables = True
900
901   if binaries_present(master_binding_instance, instances) is True:
902
903      print >>output_file, "binary"
904      print >>output_file, ""
905
906      print >>output_file, "\\ Begin binary variable block for master"
907      problem_writer._print_model_LP(master_binding_instance, output_file)
908      print >>output_file, "\\ End binary variable block for master"
909      print >>output_file, ""
910   
911      for scenario_name in instances.keys():
912         instance = instances[scenario_name]
913         print >>output_file, "\\ Begin binary variable block for scenario",scenario_name
914         problem_writer._print_model_LP(instance, output_file)
915         print >>output_file, "\\ End integer binary block for scenario",scenario_name
916         print >>output_file, ""
917
918   # wrap up.
919   print >>output_file, "end"
920
921   # clean up.
922   output_file.close()
923
924   print ""
925   print "Output file written to file=",output_filename
926
927   print ""
928   print "Done..."
929
930   end_time = time.time()
931
932   print ""
933   print "Total execution time=%8.2f seconds" %(end_time - start_time)
934   print ""   
Note: See TracBrowser for help on using the repository browser.