source: coopr.pysp/stable/2.1/coopr/pysp/ph.py @ 1952

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

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

........

r1884 | wehart | 2009-11-14 00:11:40 -0700 (Sat, 14 Nov 2009) | 2 lines


Update to current revision number.

........

r1890 | jwatson | 2009-11-14 17:08:09 -0700 (Sat, 14 Nov 2009) | 3 lines


Added option to PH to fix all discrete variables that are converged upon termination. Functionality is within the WW PH extension. Not for typical use, but useful for debugging.

........

r1895 | jwatson | 2009-11-16 10:44:48 -0700 (Mon, 16 Nov 2009) | 3 lines


Changed PySP package manager information from Bill to JPW.

........

r1911 | jwatson | 2009-11-19 10:24:07 -0700 (Thu, 19 Nov 2009) | 3 lines


Restructured PH code to allow for general breakpoint distribution strategies. Added option to runph (--breakpoint-strategy) to specify particular strategies.

........

r1912 | jwatson | 2009-11-19 13:24:04 -0700 (Thu, 19 Nov 2009) | 3 lines


Added PH breakpoint distribution strategy that uniformly concentrates pieces between the observed node-min and node-max values, as opposed to the full lb/ub range.

........

r1913 | jwatson | 2009-11-19 13:46:13 -0700 (Thu, 19 Nov 2009) | 3 lines


Added Woodruff PH breakpoint distribution strategy.

........

r1914 | jwatson | 2009-11-19 15:22:56 -0700 (Thu, 19 Nov 2009) | 3 lines


Added exponentially biased breakpoint distribution technique to PH.

........

r1922 | jwatson | 2009-11-20 15:50:27 -0700 (Fri, 20 Nov 2009) | 3 lines


Fix to make PH compliant with recent solution re-structuring.

........

r1925 | jwatson | 2009-11-20 22:42:49 -0700 (Fri, 20 Nov 2009) | 3 lines


Fixes to make PySP consistent with the variable lower/upper bound changes described in the previous commit. Fixed a few bugs in the farmer example now that bounds specified through parameters are working correctly.

........

r1928 | jwatson | 2009-11-21 11:33:10 -0700 (Sat, 21 Nov 2009) | 3 lines


Fixed an issue with construction of piecewise linear breakpoints in PH - if the average was at a bound, divide-by-zero was triggering.

........

r1929 | jwatson | 2009-11-21 11:45:34 -0700 (Sat, 21 Nov 2009) | 3 lines


Fixed another numerical issue in breakpoint computation in PH.

........

r1930 | jwatson | 2009-11-21 12:53:31 -0700 (Sat, 21 Nov 2009) | 3 lines


Adding 288 scenario forestry instance, for a serious challenge!

........

r1931 | jwatson | 2009-11-21 13:01:18 -0700 (Sat, 21 Nov 2009) | 3 lines


Added a feasible/new chile 48 instance (the old, hand-constructed instance seemed to be infeasible), and the EF for the 288 scenario instance.

........

r1932 | jwatson | 2009-11-21 13:37:45 -0700 (Sat, 21 Nov 2009) | 3 lines


Bug fix to PH when LB/UB aren't specified and you are linearizing.

........

r1933 | jwatson | 2009-11-21 14:24:35 -0700 (Sat, 21 Nov 2009) | 3 lines


Updating documentation on forestry PySP example.

........

r1934 | jwatson | 2009-11-22 15:38:31 -0700 (Sun, 22 Nov 2009) | 3 lines


Fix to profiler option in PH.

........

File size: 78.8 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2008 Sandia Corporation.
5#  This software is distributed under the BSD License.
6#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7#  the U.S. Government retains certain rights in this software.
8#  For more information, see the Coopr README.txt file.
9#  _________________________________________________________________________
10
11import sys
12import types
13from coopr.pyomo import *
14import copy
15import os.path
16import traceback
17import copy
18from coopr.opt import SolverResults,SolverStatus
19from coopr.opt.base import SolverFactory
20from coopr.opt.parallel import SolverManagerFactory
21import time
22import types
23import pickle
24from math import fabs, log, exp
25
26from scenariotree import *
27from phutils import *
28
29from pyutilib.plugin.core import ExtensionPoint
30
31from coopr.pysp.phextension import IPHExtension
32
33class ProgressiveHedging(object):
34
35   #
36   # routine to compute linearization breakpoints uniformly between the bounds and the mean.
37   #
38   
39   def compute_uniform_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
40     
41      breakpoints = []
42
43      # add the lower bound - the first breakpoint.
44      breakpoints.append(lb)
45
46      # determine the breakpoints to the left of the mean.
47      left_step = (xavg - lb) / num_breakpoints_per_side
48      current_x = lb
49      for i in range(1,num_breakpoints_per_side+1):
50         this_lb = current_x
51         this_ub = current_x+left_step
52         if (fabs(this_lb-lb) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
53            breakpoints.append(this_lb)
54         current_x += left_step           
55
56      # add the mean - it's always a breakpoint. unless!
57      # the lb or ub and the avg are the same.
58      if (fabs(lb-xavg) > self._integer_tolerance) and (fabs(ub-xavg) > self._integer_tolerance):
59         breakpoints.append(xavg)
60
61      # determine the breakpoints to the right of the mean.
62      right_step = (ub - xavg) / num_breakpoints_per_side
63      current_x = xavg
64      for i in range(1,num_breakpoints_per_side+1):
65         this_lb = current_x
66         this_ub = current_x+right_step
67         if (fabs(this_ub-xavg) > self._integer_tolerance) and (fabs(this_ub-ub) > self._integer_tolerance):         
68            breakpoints.append(this_ub)
69         current_x += right_step
70
71      # add the upper bound - the last breakpoint.
72      # the upper bound should always be different than the lower bound (I say with some
73      # hesitation - it really depends on what plugins are doing to modify the bounds dynamically).
74      breakpoints.append(ub)
75
76      return breakpoints
77
78   #
79   # routine to compute linearization breakpoints uniformly between the current node min/max bounds.
80   #   
81
82   def compute_uniform_between_nodestat_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
83
84      breakpoints = []
85
86      # add the lower bound - the first breakpoint.
87      breakpoints.append(lb)
88
89      # add the node-min - the second breakpoint. but only if it is different than the lower bound and the mean.
90      if (fabs(node_min-lb) > self._integer_tolerance) and (fabs(node_min-xavg) > self._integer_tolerance):     
91         breakpoints.append(node_min)
92
93      step = (node_max - node_min) / num_breakpoints
94      current_x = node_min
95      for i in range(1,num_breakpoints+1):
96         this_lb = current_x
97         this_ub = current_x+step
98         if (fabs(this_lb-node_min) > self._integer_tolerance) and (fabs(this_lb-node_max) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
99            breakpoints.append(this_lb)
100         current_x += step           
101
102      # add the node-max - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
103      if (fabs(node_max-ub) > self._integer_tolerance) and (fabs(node_max-xavg) > self._integer_tolerance):           
104         breakpoints.append(node_max)     
105
106      # add the upper bound - the last breakpoint.
107      breakpoints.append(ub)
108
109      # add the mean - it's always a breakpoint. unless!
110      # it happens to be equal to (within tolerance) the lower or upper bounds.
111      # sort to insert it in the correct place.
112      if (fabs(xavg - lb) > self._integer_tolerance) and (fabs(xavg - ub) > self._integer_tolerance):
113         breakpoints.append(xavg)
114      breakpoints.sort()
115
116      return breakpoints
117
118   #
119   # routine to compute linearization breakpoints using "Woodruff" relaxation of the compute_uniform_between_nodestat_breakpoints.
120   #   
121
122   def compute_uniform_between_woodruff_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints):
123
124      breakpoints = []
125
126      # add the lower bound - the first breakpoint.
127      breakpoints.append(lb)
128
129      # be either three "differences" from the mean, or else "halfway to the bound", whichever is closer to the mean.
130      left = max(xavg - 3.0 * (xavg - node_min), xavg - 0.5 * (xavg - lb))
131      right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))     
132
133      # add the left bound - the second breakpoint. but only if it is different than the lower bound and the mean.
134      if (fabs(left-lb) > self._integer_tolerance) and (fabs(left-xavg) > self._integer_tolerance):     
135         breakpoints.append(left)
136
137      step = (right - left) / num_breakpoints
138      current_x = left
139      for i in range(1,num_breakpoints+1):
140         this_lb = current_x
141         this_ub = current_x+step
142         if (fabs(this_lb-left) > self._integer_tolerance) and (fabs(this_lb-right) > self._integer_tolerance) and (fabs(this_lb-xavg) > self._integer_tolerance):
143            breakpoints.append(this_lb)
144         current_x += step           
145
146      # add the right bound - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
147      if (fabs(right-ub) > self._integer_tolerance) and (fabs(right-xavg) > self._integer_tolerance):           
148         breakpoints.append(right)     
149
150      # add the upper bound - the last breakpoint.
151      breakpoints.append(ub)
152
153      # add the mean - it's always a breakpoint.
154      # sort to insert it in the correct place.
155      breakpoints.append(xavg)
156      breakpoints.sort()
157
158      return breakpoints
159
160   #
161   # routine to compute linearization breakpoints based on an exponential distribution from the mean in each direction.
162   #   
163
164   def compute_exponential_from_mean_breakpoints(self, lb, node_min, xavg, node_max, ub, num_breakpoints_per_side):
165
166      breakpoints = []
167
168      # add the lower bound - the first breakpoint.
169      breakpoints.append(lb)
170
171      # determine the breakpoints to the left of the mean.
172      left_delta = xavg - lb
173      base = exp(log(left_delta) / num_breakpoints_per_side)
174      current_offset = base
175      for i in range(1,num_breakpoints_per_side+1):
176         current_x = xavg - current_offset
177         if (fabs(current_x-lb) > self._integer_tolerance) and (fabs(current_x-xavg) > self._integer_tolerance):
178            breakpoints.append(current_x)
179         current_offset *= base
180
181      # add the mean - it's always a breakpoint.
182      breakpoints.append(xavg)
183
184      # determine the breakpoints to the right of the mean.
185      right_delta = ub - xavg
186      base = exp(log(right_delta) / num_breakpoints_per_side)
187      current_offset = base
188      for i in range(1,num_breakpoints_per_side+1):
189         current_x = xavg + current_offset
190         if (fabs(current_x-xavg) > self._integer_tolerance) and (fabs(current_x-ub) > self._integer_tolerance):         
191            breakpoints.append(current_x)
192         current_offset *= base
193
194      # add the upper bound - the last breakpoint.
195      breakpoints.append(ub)         
196
197      return breakpoints     
198
199   #
200   # a utility intended for folks who are brave enough to script rho setting in a python file.
201   #
202
203   def setRhoAllScenarios(self, variable_value, rho_expression):
204
205      variable_name = None
206      variable_index = None
207
208      if isVariableNameIndexed(variable_value.name) is True:
209
210         variable_name, variable_index = extractVariableNameAndIndex(variable_value.name)
211
212      else:
213
214         variable_name = variable_value.name
215         variable_index = None
216     
217      new_rho_value = rho_expression()
218
219      if self._verbose is True:
220         print "Setting rho="+str(new_rho_value)+" for variable="+variable_value.name
221
222      for instance_name, instance in self._instances.items():
223
224         rho_param = getattr(instance, "PHRHO_"+variable_name)
225         rho_param[variable_index] = new_rho_value
226
227   #
228   # a utility intended for folks who are brave enough to script variable bounds setting in a python file.
229   #
230
231   def setVariableBoundsAllScenarios(self, variable_name, variable_index, lower_bound, upper_bound):
232
233      if isinstance(lower_bound, float) is False:
234         raise ValueError, "Lower bound supplied to PH method setVariableBoundsAllScenarios for variable="+variable_name+indexToString(variable_index)+" must be a constant; value supplied="+str(lower_bound)
235
236      if isinstance(upper_bound, float) is False:
237         raise ValueError, "Upper bound supplied to PH method setVariableBoundsAllScenarios for variable="+variable_name+indexToString(variable_index)+" must be a constant; value supplied="+str(upper_bound)
238
239      for instance_name, instance in self._instances.items():
240
241         variable = getattr(instance, variable_name)
242         variable[variable_index].setlb(lower_bound)
243         variable[variable_index].setub(upper_bound)
244
245   #
246   # a utility intended for folks who are brave enough to script variable bounds setting in a python file.
247   # same functionality as above, but applied to all indicies of the variable, in all scenarios.
248   #
249
250   def setVariableBoundsAllIndicesAllScenarios(self, variable_name, lower_bound, upper_bound):
251
252      if isinstance(lower_bound, float) is False:
253         raise ValueError, "Lower bound supplied to PH method setVariableBoundsAllIndiciesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(lower_bound)
254
255      if isinstance(upper_bound, float) is False:
256         raise ValueError, "Upper bound supplied to PH method setVariableBoundsAllIndicesAllScenarios for variable="+variable_name+" must be a constant; value supplied="+str(upper_bound)
257
258      for instance_name, instance in self._instances.items():
259
260         variable = getattr(instance, variable_name)
261         for index in variable:
262            variable[index].setlb(lower_bound)
263            variable[index].setub(upper_bound)
264
265   #
266   # checkpoint the current PH state via pickle'ing. the input iteration count
267   # simply serves as a tag to create the output file name. everything with the
268   # exception of the _ph_plugins, _solver_manager, and _solver attributes are
269   # pickled. currently, plugins fail in the pickle process, which is fine as
270   # JPW doesn't think you want to pickle plugins (particularly the solver and
271   # solver manager) anyway. For example, you might want to change those later,
272   # after restoration - and the PH state is independent of how scenario
273   # sub-problems are solved.
274   #
275
276   def checkpoint(self, iteration_count):
277
278      checkpoint_filename = "checkpoint."+str(iteration_count)
279
280      tmp_ph_plugins = self._ph_plugins
281      tmp_solver_manager = self._solver_manager
282      tmp_solver = self._solver
283
284      self._ph_plugins = None
285      self._solver_manager = None
286      self._solver = None
287     
288      checkpoint_file = open(checkpoint_filename, "w")
289      pickle.dump(self,checkpoint_file)
290      checkpoint_file.close()
291
292      self._ph_plugins = tmp_ph_plugins
293      self._solver_manager = tmp_solver_manager
294      self._solver = tmp_solver
295
296      print "Checkpoint written to file="+checkpoint_filename
297   
298   #
299   # a simple utility to count the number of continuous and discrete variables in a set of instances.
300   # unused variables are ignored, and counts include all active indices. returns a pair - num-discrete,
301   # num-continuous.
302   #
303
304   def compute_variable_counts(self):
305
306      num_continuous_vars = 0
307      num_discrete_vars = 0
308     
309      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
310
311         for tree_node in stage._tree_nodes:
312
313            for (variable, index_template, variable_indices) in stage._variables:
314
315               for index in variable_indices:
316
317                  is_used = True # until proven otherwise                     
318                  for scenario in tree_node._scenarios:
319                     instance = self._instances[scenario._name]
320                     if getattr(instance,variable.name)[index].status == VarStatus.unused:
321                        is_used = False
322
323                  if is_used is True:                       
324
325                     if isinstance(variable.domain, IntegerSet) or isinstance(variable.domain, BooleanSet):
326                        num_discrete_vars = num_discrete_vars + 1
327                     else:
328                        num_continuous_vars = num_continuous_vars + 1
329
330      return (num_discrete_vars, num_continuous_vars)
331
332   #
333   # ditto above, but count the number of fixed discrete and continuous variables.
334   # important: once a variable (value) is fixed, it is flagged as unused in the
335   # course of presolve - because it is no longer referenced. this makes sense,
336   # of course; it's just something to watch for. this is an obvious assumption
337   # that we won't be fixing unused variables, which should not be an issue.
338   #
339
340   def compute_fixed_variable_counts(self):
341
342      num_fixed_continuous_vars = 0
343      num_fixed_discrete_vars = 0
344     
345      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
346         
347         for tree_node in stage._tree_nodes:
348
349            for (variable, index_template, variable_indices) in stage._variables:
350
351               for index in variable_indices:
352
353                  # implicit assumption is that if a variable value is fixed in one
354                  # scenario, it is fixed in all scenarios.
355
356                  is_fixed = False # until proven otherwise
357                  for scenario in tree_node._scenarios:
358                     instance = self._instances[scenario._name]
359                     var_value = getattr(instance,variable.name)[index]
360                     if var_value.fixed is True:
361                        is_fixed = True
362
363                  if is_fixed is True:
364
365                     if isinstance(variable.domain, IntegerSet) or isinstance(variable.domain, BooleanSet):
366                        num_fixed_discrete_vars = num_fixed_discrete_vars + 1
367                     else:
368                        num_fixed_continuous_vars = num_fixed_continuous_vars + 1                           
369
370      return (num_fixed_discrete_vars, num_fixed_continuous_vars)
371
372   # a utility to create piece-wise linear constraint expressions for a given variable, for
373   # use in constructing the augmented (penalized) PH objective. lb and ub are the bounds
374   # on this piece, variable is the actual instance variable, and average is the instance
375   # parameter specifying the average of this variable across instances sharing passing
376   # through a common tree node. lb and ub are floats.
377   def _create_piecewise_constraint_expression(self, lb, ub, instance_variable, variable_average, quad_variable):
378
379      penalty_at_lb = (lb - variable_average()) * (lb - variable_average())
380      penalty_at_ub = (ub - variable_average()) * (ub - variable_average())
381      slope = (penalty_at_ub - penalty_at_lb) / (ub - lb)
382      intercept = penalty_at_lb - slope * lb
383      expression = (0.0, quad_variable - slope * instance_variable - intercept, None)
384
385      return expression
386
387   # when the quadratic penalty terms are approximated via piecewise linear segments,
388   # we end up (necessarily) "littering" the scenario instances with extra constraints.
389   # these need to and should be cleaned up after PH, for purposes of post-PH manipulation,
390   # e.g., writing the extensive form.
391   def _cleanup_scenario_instances(self):
392
393      for instance_name, instance in self._instances.items():
394
395         for constraint_name in self._instance_augmented_attributes[instance_name]:
396
397            instance._clear_attribute(constraint_name)
398
399         # if you don't pre-solve, the name collections won't be updated.
400         instance.presolve()
401
402   # create PH weight and xbar vectors, on a per-scenario basis, for each variable that is not in the
403   # final stage, i.e., for all variables that are being blended by PH. the parameters are created
404   # in the space of each scenario instance, so that they can be directly and automatically
405   # incorporated into the (appropriately modified) objective function.
406
407   def _create_ph_scenario_parameters(self):
408
409      for (instance_name, instance) in self._instances.items():
410
411         # first, gather all unique variables referenced in any stage
412         # other than the last, independent of specific indices. this
413         # "gather" step is currently required because we're being lazy
414         # in terms of index management in the scenario tree - which
415         # should really be done in terms of sets of indices.
416         # NOTE: technically, the "instance variables" aren't really references
417         # to the variable in the instance - instead, the reference model. this
418         # isn't an issue now, but it could easily become one (esp. in avoiding deep copies).
419         instance_variables = {}
420         
421         for stage in self._scenario_tree._stages[:-1]:
422           
423            for (reference_variable, index_template, reference_indices) in stage._variables:
424                 
425               if reference_variable.name not in instance_variables.keys():
426                     
427                  instance_variables[reference_variable.name] = reference_variable
428
429         # for each blended variable, create a corresponding ph weight and average parameter in the instance.
430         # this is a bit wasteful, in terms of indices that might appear in the last stage, but that is minor
431         # in the grand scheme of things.
432
433         for (variable_name, reference_variable) in instance_variables.items():
434
435            # PH WEIGHT
436
437            new_w_index = reference_variable._index
438            new_w_parameter_name = "PHWEIGHT_"+reference_variable.name
439            new_w_parameter = Param(new_w_index,name=new_w_parameter_name)
440            setattr(instance,new_w_parameter_name,new_w_parameter)
441
442            # if you don't explicitly assign values to each index, the entry isn't created - instead, when you reference
443            # the parameter that hasn't been explicitly assigned, you just get the default value as a constant. I'm not
444            # sure if this has to do with the model output, or the function of the model, but I'm doing this to avoid the
445            # issue in any case for now.
446            for index in new_w_index:
447               new_w_parameter[index] = 0.0
448
449            # PH AVG
450
451            new_avg_index = reference_variable._index
452            new_avg_parameter_name = "PHAVG_"+reference_variable.name
453            new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
454            setattr(instance,new_avg_parameter_name,new_avg_parameter)
455
456            for index in new_avg_index:
457               new_avg_parameter[index] = 0.0
458
459            # PH RHO
460
461            new_rho_index = reference_variable._index
462            new_rho_parameter_name = "PHRHO_"+reference_variable.name
463            new_rho_parameter = Param(new_rho_index,name=new_rho_parameter_name)
464            setattr(instance,new_rho_parameter_name,new_rho_parameter)
465
466            for index in new_rho_index:
467               new_rho_parameter[index] = self._rho
468
469            # PH LINEARIZED PENALTY TERM
470
471            if self._linearize_nonbinary_penalty_terms > 0:
472
473               new_penalty_term_variable_index = reference_variable._index
474               new_penalty_term_variable_name = "PHQUADPENALTY_"+reference_variable.name
475               # this is a quadratic penalty term - the lower bound is 0!
476               new_penalty_term_variable = Var(new_penalty_term_variable_index, name=new_penalty_term_variable_name, bounds=(0.0,None))
477               new_penalty_term_variable._construct(instance,None)                                             
478               setattr(instance, new_penalty_term_variable_name, new_penalty_term_variable)
479               self._instance_augmented_attributes[instance_name].append(new_penalty_term_variable_name)
480
481            # BINARY INDICATOR PARAMETER FOR WHETHER SPECIFIC VARIABLES ARE BLENDED. FOR ADVANCED USERS ONLY.
482
483            # also controls whether weight updates proceed at any iteration.
484
485            new_blend_index = reference_variable._index
486            new_blend_parameter_name = "PHBLEND_"+reference_variable.name
487            new_blend_parameter = Param(new_blend_index, name=new_blend_parameter_name, within=Binary)
488            setattr(instance,new_blend_parameter_name,new_blend_parameter)
489
490            for index in new_blend_index:
491               new_blend_parameter[index] = 1
492     
493
494   """ Constructor
495       Arguments:
496          max_iterations        the maximum number of iterations to run PH (>= 0). defaults to 0.
497          rho                   the global rho value (> 0). defaults to 0.
498          rho_setter            an optional name of a python file used to set particular variable rho values.
499          solver                the solver type that PH uses to solve scenario sub-problems. defaults to "cplex".
500          solver_manager        the solver manager type that coordinates scenario sub-problem solves. defaults to "serial".
501          keep_solver_files     do I keep intermediate solver files around (for debugging)? defaults to False.
502          output_solver_log     do I dump the solver log (as it is being generated) to the screen? defaults to False.
503          output_solver_results do I output (for debugging) the detailed solver results, including solutions, for scenario solves? defaults to False.
504          verbose               does the PH object stream debug/status output? defaults to False.
505          output_times          do I output timing statistics? defaults to False (e.g., useful in the case where you want to regression test against baseline output).
506          checkpoint_interval   how many iterations between writing a checkpoint file containing the entire PH state? defaults to 0, indicating never.
507
508   """
509   def __init__(self, *args, **kwds):
510
511      # PH configuration parameters
512      self._rho = 0.0 # a default, global values for rhos.
513      self._rho_setter = None # filename for the modeler to set rho on a per-variable or per-scenario basis.
514      self._bounds_setter = None # filename for the modeler to set rho on a per-variable basis, after all scenarios are available.
515      self._max_iterations = 0
516
517      # PH reporting parameters
518      self._verbose = False # do I flood the screen with status output?
519      self._output_continuous_variable_stats = True # when in verbose mode, do I output weights/averages for continuous variables?
520      self._output_solver_results = False
521      self._output_times = False
522
523      # PH run-time variables
524      self._current_iteration = 0 # the 'k'
525      self._xbar = {} # current per-variable averages. maps (node_id, variable_name) -> value
526      self._initialized = False # am I ready to call "solve"? Set to True by the initialize() method.
527
528      # PH solver information / objects.
529      self._solver_type = "cplex"
530      self._solver_manager_type = "serial" # serial or pyro are the options currently available
531     
532      self._solver = None # will eventually be unnecessary once Bill eliminates the need for a solver in the solver manager constructor.
533      self._solver_manager = None 
534     
535      self._keep_solver_files = False
536      self._output_solver_log = False
537
538      # PH convergence computer/updater.
539      self._converger = None
540
541      # PH history
542      self._solutions = {}
543
544      # the checkpoint interval - expensive operation, but worth it for big models.
545      # 0 indicates don't checkpoint.
546      self._checkpoint_interval = 0
547
548      # all information related to the scenario tree (implicit and explicit).
549      self._model = None # not instantiated
550      self._model_instance = None # instantiated
551
552      self._scenario_tree = None
553     
554      self._scenario_data_directory = "" # this the prefix for all scenario data
555      self._instances = {} # maps scenario name to the corresponding model instance
556      self._instance_augmented_attributes = {} # maps scenario name to a list of the constraints added (e.g., for piecewise linear approximation) to the instance by PH.
557
558      # for various reasons (mainly hacks at this point), it's good to know whether we're minimizing or maximizing.
559      self._is_minimizing = None
560
561      # global handle to ph extension plugins
562      self._ph_plugins = ExtensionPoint(IPHExtension)
563
564      # PH timing statistics - relative to last invocation.
565      self._init_start_time = None # for initialization() method
566      self._init_end_time = None
567      self._solve_start_time = None # for solve() method
568      self._solve_end_time = None
569      self._cumulative_solve_time = None # seconds, over course of solve()
570      self._cumulative_xbar_time = None # seconds, over course of update_xbars()
571      self._cumulative_weight_time = None # seconds, over course of update_weights()
572
573      # do I disable warm-start for scenario sub-problem solves during PH iterations >= 1?
574      self._disable_warmstarts = False
575
576      # do I drop proximal (quadratic penalty) terms from the weighted objective functions?
577      self._drop_proximal_terms = False
578
579      # do I linearize the quadratic penalty term for continuous variables via a
580      # piecewise linear approximation? the default should always be 0 (off), as the
581      # user should be aware when they are forcing an approximation.
582      self._linearize_nonbinary_penalty_terms = 0
583
584      # the breakpoint distribution strategy employed when linearizing. 0 implies uniform
585      # distribution between the variable lower and upper bounds.
586      self._breakpoint_strategy = 0
587
588      # do I retain quadratic objective terms associated with binary variables? in general,
589      # there is no good reason to not linearize, but just in case, I introduced the option.
590      self._retain_quadratic_binary_terms = False
591
592      # PH default tolerances - for use in fixing and testing equality across scenarios,
593      # and other stuff.
594      self._integer_tolerance = 0.00001
595
596      # process the keyword options
597      for key in kwds.keys():
598         if key == "max_iterations":
599            self._max_iterations = kwds[key]
600         elif key == "rho":
601            self._rho = kwds[key]
602         elif key == "rho_setter":
603            self._rho_setter = kwds[key]
604         elif key == "bounds_setter":
605            self._bounds_setter = kwds[key]                       
606         elif key == "solver":
607            self._solver_type = kwds[key]
608         elif key == "solver_manager":
609            self._solver_manager_type = kwds[key]           
610         elif key == "keep_solver_files":
611            self._keep_solver_files = kwds[key]
612         elif key == "output_solver_results":
613            self._output_solver_results = kwds[key]
614         elif key == "output_solver_log":
615            self._output_solver_log = kwds[key]
616         elif key == "verbose":
617            self._verbose = kwds[key]
618         elif key == "output_times":
619            self._output_times = kwds[key]
620         elif key == "disable_warmstarts":
621            self._disable_warmstarts = kwds[key]
622         elif key == "drop_proximal_terms":
623            self._drop_proximal_terms = kwds[key]
624         elif key == "retain_quadratic_binary_terms":
625            self._retain_quadratic_binary_terms = kwds[key]
626         elif key == "linearize_nonbinary_penalty_terms":
627            self._linearize_nonbinary_penalty_terms = kwds[key]
628         elif key == "breakpoint_strategy":
629            self._breakpoint_strategy = kwds[key]
630         elif key == "checkpoint_interval":
631            self._checkpoint_interval = kwds[key]           
632         else:
633            print "Unknown option=" + key + " specified in call to PH constructor"
634
635      # validate all "atomic" options (those that can be validated independently)
636      if self._max_iterations < 0:
637         raise ValueError, "Maximum number of PH iterations must be non-negative; value specified=" + `self._max_iterations`
638      if self._rho <= 0.0:
639         raise ValueError, "Value of rho paramter in PH must be non-zero positive; value specified=" + `self._rho`
640
641      # validate the linearization (number of pieces) and breakpoint distribution parameters.
642      if self._linearize_nonbinary_penalty_terms < 0:
643         raise ValueError, "Value of linearization parameter for nonbinary penalty terms must be non-negative; value specified=" + `self._linearize_nonbinary_penalty_terms`         
644      if self._breakpoint_strategy < 0:
645         raise ValueError, "Value of the breakpoint distribution strategy parameter must be non-negative; value specified=" + str(self._breakpoint_strategy)
646      if self._breakpoint_strategy > 3:
647         raise ValueError, "Unknown breakpoint distribution strategy specified - valid values are between 0 and 2, inclusive; value specified=" + str(self._breakpoint_strategy)
648   
649      # validate rho setter file if specified.
650      if self._rho_setter is not None:
651         if os.path.exists(self._rho_setter) is False:
652            raise ValueError, "The rho setter script file="+self._rho_setter+" does not exist"
653
654      # validate bounds setter file if specified.
655      if self._bounds_setter is not None:
656         if os.path.exists(self._bounds_setter) is False:
657            raise ValueError, "The bounds setter script file="+self._bounds_setter+" does not exist"     
658
659      # validate the checkpoint interval.
660      if self._checkpoint_interval < 0:
661         raise ValueError, "A negative checkpoint interval with value="+str(self._checkpoint_interval)+" was specified in call to PH constructor"
662
663      # construct the sub-problem solver.
664      if self._verbose is True:
665         print "Constructing solver type="+self._solver_type         
666      self._solver = SolverFactory(self._solver_type)
667      if self._solver == None:
668         raise ValueError, "Unknown solver type=" + self._solver_type + " specified in call to PH constructor"
669      if self._keep_solver_files is True:
670         self._solver.keepFiles = True
671
672      # construct the solver manager.
673      if self._verbose is True:
674         print "Constructing solver manager of type="+self._solver_manager_type
675      self._solver_manager = SolverManagerFactory(self._solver_manager_type)
676      if self._solver_manager is None:
677         raise ValueError, "Failed to create solver manager of type="+self._solver_manager_type+" specified in call to PH constructor"
678
679      # a set of all valid PH iteration indicies is generally useful for plug-ins, so create it here.
680      self._iteration_index_set = Set(name="PHIterations")
681      for i in range(0,self._max_iterations + 1):
682         self._iteration_index_set.add(i)     
683
684      # spit out parameterization if verbosity is enabled
685      if self._verbose is True:
686         print "PH solver configuration: "
687         print "   Max iterations=" + `self._max_iterations`
688         print "   Default global rho=" + `self._rho`
689         if self._rho_setter is not None:
690            print "   Rho initialization file=" + self._rho_setter
691         if self._bounds_setter is not None:
692            print "   Variable bounds initialization file=" + self._bounds_setter           
693         print "   Sub-problem solver type=" + `self._solver_type`
694         print "   Solver manager type=" + `self._solver_manager_type`
695         print "   Keep solver files? " + str(self._keep_solver_files)
696         print "   Output solver results? " + str(self._output_solver_results)
697         print "   Output solver log? " + str(self._output_solver_log)
698         print "   Output times? " + str(self._output_times)
699         print "   Checkpoint interval="+str(self._checkpoint_interval)
700
701   """ Initialize PH with model and scenario data, in preparation for solve().
702       Constructs and reads instances.
703   """
704   def initialize(self, scenario_data_directory_name=".", model=None, model_instance=None, scenario_tree=None, converger=None):
705
706      self._init_start_time = time.time()
707
708      if self._verbose is True:
709         print "Initializing PH"
710         print "   Scenario data directory=" + scenario_data_directory_name
711
712      if not os.path.exists(scenario_data_directory_name):
713         raise ValueError, "Scenario data directory=" + scenario_data_directory_name + " either does not exist or cannot be read"
714
715      self._scenario_data_directory_name = scenario_data_directory_name
716
717      # IMPT: The input model should be an *instance*, as it is very useful (critical!) to know
718      #       the dimensions of sets, be able to store suffixes on variable values, etc.
719      if model is None:
720         raise ValueError, "A model must be supplied to the PH initialize() method"
721
722      if scenario_tree is None:
723         raise ValueError, "A scenario tree must be supplied to the PH initialize() method"
724
725      if converger is None:
726         raise ValueError, "A convergence computer must be supplied to the PH initialize() method"
727
728      self._model = model
729      self._model_instance = model_instance
730      self._scenario_tree = scenario_tree
731      self._converger = converger
732
733      model_objective = model._component[Objective]
734      self._is_minimizing = (model_objective[ model_objective.keys()[0] ].sense == minimize)
735
736      self._converger.reset()
737
738      # construct the instances for each scenario
739      if self._verbose is True:
740         if self._scenario_tree._scenario_based_data == 1:
741            print "Scenario-based instance initialization enabled"
742         else:
743            print "Node-based instance initialization enabled"
744         
745      for scenario in self._scenario_tree._scenarios:
746
747         scenario_instance = None
748
749         if self._verbose is True:
750            print "Creating instance for scenario=" + scenario._name
751
752         try:
753            if self._scenario_tree._scenario_based_data == 1:
754               scenario_data_filename = self._scenario_data_directory_name + os.sep + scenario._name + ".dat"
755               if self._verbose is True:
756                  print "Data for scenario=" + scenario._name + " loads from file=" + scenario_data_filename
757               scenario_instance = (self._model).create(scenario_data_filename)
758            else:
759               scenario_instance = self._model.clone()
760               scenario_data = ModelData()
761               current_node = scenario._leaf_node
762               while current_node is not None:
763                  node_data_filename = self._scenario_data_directory_name + os.sep + current_node._name + ".dat"
764                  if self._verbose is True:
765                     print "Node data for scenario=" + scenario._name + " partially loading from file=" + node_data_filename
766                  scenario_data.add_data_file(node_data_filename)
767                  current_node = current_node._parent
768               scenario_data.read(model=scenario_instance)
769               scenario_instance._load_model_data(scenario_data)
770               scenario_instance.presolve()
771         except:
772            print "Encountered exception in model instance creation - traceback:"
773            traceback.print_exc()
774            raise RuntimeError, "Failed to create model instance for scenario=" + scenario._name
775
776         self._instances[scenario._name] = scenario_instance
777         self._instances[scenario._name].name = scenario._name
778         self._instance_augmented_attributes[scenario._name] = []
779
780      # create ph-specific parameters (weights, xbar, etc.) for each instance.
781
782      if self._verbose is True:
783         print "Creating weight, average, and rho parameter vectors for scenario instances"
784
785      self._create_ph_scenario_parameters()
786           
787      # if specified, run the user script to initialize variable rhos at their whim.
788      if self._rho_setter is not None:
789         print "Executing user rho set script from filename=", self._rho_setter
790         execfile(self._rho_setter)
791
792      # with the instances created, run the user script to initialize variable bounds.
793      if self._bounds_setter is not None:
794         print "Executing user variable bounds set script from filename=", self._bounds_setter
795         execfile(self._bounds_setter)
796
797      # create parameters to store variable statistics (of general utility) at each node in the scenario tree.
798
799      if self._verbose is True:
800         print "Creating variable statistic (min/avg/max) parameter vectors for scenario tree nodes"
801
802      for stage in self._scenario_tree._stages[:-1]:
803
804         # first, gather all unique variables referenced in this stage
805         # this "gather" step is currently required because we're being lazy
806         # in terms of index management in the scenario tree - which
807         # should really be done in terms of sets of indices.
808
809         stage_variables = {}
810         for (reference_variable, index_template, reference_index) in stage._variables:
811            if reference_variable.name not in stage_variables.keys():
812               stage_variables[reference_variable.name] = reference_variable
813
814         # next, create min/avg/max parameters for each variable in the corresponding tree node.
815         # NOTE: the parameter names below could really be empty, as they are never referenced
816         #       explicitly.
817         for (variable_name, reference_variable) in stage_variables.items():
818            for tree_node in stage._tree_nodes:
819
820               new_min_index = reference_variable._index
821               new_min_parameter_name = "NODEMIN_"+reference_variable.name
822               new_min_parameter = Param(new_min_index,name=new_min_parameter_name)
823               for index in new_min_index:
824                  new_min_parameter[index] = 0.0
825               tree_node._minimums[reference_variable.name] = new_min_parameter                                   
826
827               new_avg_index = reference_variable._index
828               new_avg_parameter_name = "NODEAVG_"+reference_variable.name
829               new_avg_parameter = Param(new_avg_index,name=new_avg_parameter_name)
830               for index in new_avg_index:
831                  new_avg_parameter[index] = 0.0
832               tree_node._averages[reference_variable.name] = new_avg_parameter
833
834               new_max_index = reference_variable._index
835               new_max_parameter_name = "NODEMAX_"+reference_variable.name
836               new_max_parameter = Param(new_max_index,name=new_max_parameter_name)
837               for index in new_max_index:
838                  new_max_parameter[index] = 0.0
839               tree_node._maximums[reference_variable.name] = new_max_parameter
840
841      # the objective functions are modified throughout the course of PH iterations.
842      # save the original, as a baseline to modify in subsequent iterations. reserve
843      # the original objectives, for subsequent modification.
844      self._original_objective_expression = {}
845      for instance_name, instance in self._instances.items():
846         objective_name = instance._component[Objective].keys()[0]
847         self._original_objective_expression[instance_name] = instance._component[Objective][objective_name]._data[None].expr
848
849      # cache the number of discrete and continuous variables in the master instance. this value
850      # is of general use, e.g., in the converger classes and in plugins.
851      (self._total_discrete_vars,self._total_continuous_vars) = self.compute_variable_counts()
852      if self._verbose is True:
853         print "Total number of discrete instance variables="+str(self._total_discrete_vars)
854         print "Total number of continuous instance variables="+str(self._total_continuous_vars)
855
856      # track the total number of fixed variables of each category at the end of each PH iteration.
857      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()
858
859      # indicate that we're ready to run.
860      self._initialized = True
861
862      if self._verbose is True:
863         print "PH successfully created model instances for all scenarios"
864
865      self._init_end_time = time.time()
866
867      if self._verbose is True:
868         print "PH is successfully initialized"
869         if self._output_times is True:
870            print "Initialization time=%8.2f seconds" % (self._init_end_time - self._init_start_time)
871
872      # let plugins know if they care.
873      if self._verbose is True:
874         print "Initializing PH plugins"
875      for plugin in self._ph_plugins:
876         plugin.post_ph_initialization(self)
877      if self._verbose is True:
878         print "PH plugin initialization complete"
879
880   """ Perform the non-weighted scenario solves and form the initial w and xbars.
881   """   
882   def iteration_0_solve(self):
883
884      if self._verbose is True:
885         print "------------------------------------------------"
886         print "Starting PH iteration 0 solves"
887
888      self._current_iteration = 0
889
890      solve_start_time = time.time()
891
892      # STEP 1: queue up the solves for all scenario sub-problems and
893      #         grab all the action handles for the subsequent barrier sync.
894
895      action_handles = []
896      action_handle_instance_map = {}
897
898      for scenario in self._scenario_tree._scenarios:
899         
900         instance = self._instances[scenario._name]
901         
902         if self._verbose is True:
903            print "Queuing solve for scenario=" + scenario._name
904
905         # IMPT: You have to re-presolve if approximating continuous variable penalty terms with a
906         #       piecewise linear function. otherwise, the newly introduced variables won't be flagged
907         #       as unused (as is correct for iteration 0), and the output writer will crater.
908         if self._linearize_nonbinary_penalty_terms > 0:
909           instance.presolve()           
910           
911         # there's nothing to warm-start from in iteration 0, so don't include the keyword in the solve call.
912         # the reason you don't want to include it is that some solvers don't know how to handle the keyword
913         # at all (despite it being false). you might want to solve iteration 0 solves using some other solver.
914
915         new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)
916         action_handle_instance_map[scenario._name] = new_action_handle
917
918         action_handles.append(new_action_handle)
919
920      # STEP 2: barrier sync for all scenario sub-problem solves.
921      if self._verbose is True:
922         print "Waiting for scenario sub-problem solves"
923      self._solver_manager.wait_all(action_handles)
924      if self._verbose is True:     
925         print "Scenario sub-problem solves completed"     
926
927      solve_end_time = time.time()
928      self._cumulative_solve_time += (solve_end_time - solve_start_time)
929
930      if self._output_times is True:
931         print "Aggregate sub-problem solve time=%8.2f" % (solve_end_time - solve_start_time)
932
933      # STEP 3: Load the results!
934      for scenario_name, action_handle in action_handle_instance_map.items():
935
936         if self._verbose is True:         
937            print "Successfully processed results for scenario="+scenario_name
938
939         instance = self._instances[scenario_name]
940         results = self._solver_manager.get_results(action_handle)
941
942         if len(results.solution) == 0:
943            results.write(num=1)
944            raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
945
946         if self._output_solver_results is True:
947            print "Results for scenario=",scenario_name
948            results.write(num=1)           
949
950         instance.load(results)
951
952         if self._verbose is True:                 
953            print "Successfully loaded solution for scenario="+scenario_name
954
955      if self._verbose is True:
956         print "Successfully completed PH iteration 0 solves - solution statistics:"
957         print "         Scenario              Objective                  Value"
958         for scenario in self._scenario_tree._scenarios:
959            instance = self._instances[scenario._name]
960            for objective_name in instance._component[Objective]:
961               objective = instance._component[Objective][objective_name]
962               print "%20s       %15s     %14.4f" % (scenario._name, objective.name, objective._data[None].expr())
963         print "------------------------------------------------"
964
965   #
966   # recompute the averages, minimum, and maximum statistics for all variables to be blended by PH, i.e.,
967   # not appearing in the final stage. technically speaking, the min/max aren't required by PH, but they
968   # are used often enough to warrant their computation and it's basically free if you're computing the
969   # average.
970   #
971   def update_variable_statistics(self):
972
973      start_time = time.time()
974     
975      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage
976         
977         for tree_node in stage._tree_nodes:
978               
979            for (variable, index_template, variable_indices) in stage._variables:
980                 
981               variable_name = variable.name
982                     
983               avg_parameter_name = "PHAVG_"+variable_name
984                 
985               for index in variable_indices:
986                  min = float("inf")
987                  avg = 0.0
988                  max = float("-inf")
989                  node_probability = 0.0
990                     
991                  is_used = True # until proven otherwise                     
992                  for scenario in tree_node._scenarios:
993                       
994                     instance = self._instances[scenario._name]
995                       
996                     if getattr(instance,variable_name)[index].status == VarStatus.unused:
997                        is_used = False
998                     else:                       
999                        node_probability += scenario._probability
1000                        var_value = getattr(instance, variable.name)[index].value
1001                        if var_value < min:
1002                           min = var_value
1003                        avg += (scenario._probability * var_value)
1004                        if var_value > max:
1005                           max = var_value
1006
1007                  if is_used is True:
1008                     tree_node._minimums[variable.name][index] = min
1009                     tree_node._averages[variable.name][index] = avg / node_probability
1010                     tree_node._maximums[variable.name][index] = max                       
1011
1012                     # distribute the newly computed average to the xbar variable in
1013                     # each instance/scenario associated with this node. only do this
1014                     # if the variable is used!
1015                     for scenario in tree_node._scenarios:
1016                        instance = self._instances[scenario._name]
1017                        avg_parameter = getattr(instance, avg_parameter_name)
1018                        avg_parameter[index] = avg / node_probability
1019
1020      end_time = time.time()
1021      self._cumulative_xbar_time += (end_time - start_time)
1022
1023   def update_weights(self):
1024
1025      # because the weight updates rely on the xbars, and the xbars are node-based,
1026      # I'm looping over the tree nodes and pushing weights into the corresponding scenarios.
1027      start_time = time.time()     
1028
1029      for stage in self._scenario_tree._stages[:-1]: # no blending over the final stage, so no weights to worry about.
1030         
1031         for tree_node in stage._tree_nodes:
1032
1033            for (variable, index_template, variable_indices) in stage._variables:
1034
1035               variable_name = variable.name
1036               blend_parameter_name = "PHBLEND_"+variable_name
1037               weight_parameter_name = "PHWEIGHT_"+variable_name
1038               rho_parameter_name = "PHRHO_"+variable_name
1039
1040               for index in variable_indices:
1041
1042                  tree_node_average = tree_node._averages[variable.name][index]()
1043
1044                  for scenario in tree_node._scenarios:
1045
1046                     instance = self._instances[scenario._name]
1047
1048                     if getattr(instance,variable.name)[index].status != VarStatus.unused:
1049
1050                        # we are currently not updating weights if blending is disabled for a variable.
1051                        # this is done on the premise that unless you are actively trying to move
1052                        # the variable toward the mean, the weights will blow up and be huge by the
1053                        # time that blending is activated.
1054                        variable_blend_indicator = getattr(instance, blend_parameter_name)[index]()                             
1055
1056                        # get the weight and rho parameters for this variable/index combination.
1057                        rho_value = getattr(instance, rho_parameter_name)[index]()
1058                        current_variable_weight = getattr(instance, weight_parameter_name)[index]()
1059                                               
1060                        # if I'm maximizing, invert value prior to adding (hack to implement negatives).
1061                        # probably fixed in Pyomo at this point - I just haven't checked in a long while.
1062                        if self._is_minimizing is False:
1063                           current_variable_weight = (-current_variable_weight)
1064                        current_variable_value = getattr(instance,variable.name)[index]()
1065                        new_variable_weight = current_variable_weight + variable_blend_indicator * rho_value * (current_variable_value - tree_node_average)
1066                        # I have the correct updated value, so now invert if maximizing.
1067                        if self._is_minimizing is False:
1068                           new_variable_weight = (-new_variable_weight)
1069                        getattr(instance, weight_parameter_name)[index].value = new_variable_weight
1070
1071      # we shouldn't have to re-simplify the expression, as we aren't adding any constant-variable terms - just modifying parameters.
1072
1073      end_time = time.time()
1074      self._cumulative_weight_time += (end_time - start_time)
1075
1076   def form_iteration_k_objectives(self):
1077
1078      # for each blended variable (i.e., those not appearing in the final stage),
1079      # add the linear and quadratic penalty terms to the objective.
1080      for instance_name, instance in self._instances.items():
1081         
1082         objective_name = instance._component[Objective].keys()[0]         
1083         objective = instance._component[Objective][objective_name]
1084         # clone the objective, because we're going to augment (repeatedly) the original objective.
1085         objective_expression = self._original_objective_expression[instance_name].clone() 
1086         # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
1087         quad_expression = 0.0
1088         
1089         for stage in self._scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs
1090
1091            variable_tree_node = None
1092            for node in stage._tree_nodes:
1093               for scenario in node._scenarios:
1094                  if scenario._name == instance_name:
1095                     variable_tree_node = node
1096                     break
1097
1098            for (reference_variable, index_template, variable_indices) in stage._variables:
1099
1100               variable_name = reference_variable.name
1101               variable_type = reference_variable.domain
1102
1103               w_parameter_name = "PHWEIGHT_"+variable_name
1104               w_parameter = instance._component[param._ParamBase][w_parameter_name]
1105                 
1106               average_parameter_name = "PHAVG_"+variable_name
1107               average_parameter = instance._component[param._ParamBase][average_parameter_name]
1108
1109               rho_parameter_name = "PHRHO_"+variable_name
1110               rho_parameter = instance._component[param._ParamBase][rho_parameter_name]
1111
1112               blend_parameter_name = "PHBLEND_"+variable_name
1113               blend_parameter = instance._component[param._ParamBase][blend_parameter_name]
1114
1115               node_min_parameter = variable_tree_node._minimums[variable_name]
1116               node_max_parameter = variable_tree_node._maximums[variable_name]               
1117
1118               quad_penalty_term_variable = None
1119               if self._linearize_nonbinary_penalty_terms > 0:
1120                  quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name
1121                  quad_penalty_term_variable = instance._component[var._VarBase][quad_penalty_term_variable_name]
1122
1123               instance_variable = instance._component[var._VarBase][variable_name]
1124
1125               for index in variable_indices:
1126
1127                  if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
1128
1129                     # add the linear (w-weighted) term is a consistent fashion, independent of variable type.
1130                     # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown.
1131                     objective_expression += (w_parameter[index] * instance_variable[index])
1132
1133                     # there are some cases in which a user may want to avoid the proximal term completely.
1134                     # it's of course only a good idea when there are at least bounds (both lb and ub) on
1135                     # the variables to be blended.
1136                     if self._drop_proximal_terms is False:
1137
1138                        # deal with binaries
1139                        if isinstance(variable_type, BooleanSet) is True:
1140
1141                           if self._retain_quadratic_binary_terms is False:
1142                              # this rather ugly form of the linearized quadratic expression term is required
1143                              # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing
1144                              # over the sum.
1145                              new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
1146                                         (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
1147                                         (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                             
1148                              if objective.sense is minimize:
1149                                 objective_expression += new_term
1150                              else:
1151                                 objective_expression -= new_term                                 
1152                           else:
1153                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
1154
1155                        # deal with everything else
1156                        else:
1157
1158                           if self._linearize_nonbinary_penalty_terms > 0:
1159
1160                              # the variables are easy - just a single entry.
1161                              if objective.sense is minimize:
1162                                 objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
1163                              else:
1164                                 objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
1165                                 
1166                              # the constraints are somewhat nastier.
1167
1168                              # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE"
1169                              xavg = average_parameter[index]
1170                              x = instance_variable[index]
1171
1172                              lb = None
1173                              ub = None
1174
1175                              if x.lb is None:
1176                                 raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
1177                              else:
1178                                 lb = x.lb()
1179
1180                              if x.ub is None:
1181                                 raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
1182                              else:
1183                                 ub = x.ub()
1184
1185                              node_min = node_min_parameter[index]()
1186                              node_max = node_max_parameter[index]()
1187
1188                              # compute the breakpoint sequence according to the specified strategy.
1189                              breakpoints = []
1190                              if self._breakpoint_strategy == 0:
1191                                 breakpoints = self.compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1192                              elif self._breakpoint_strategy == 1:
1193                                 breakpoints = self.compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1194                              elif self._breakpoint_strategy == 2:
1195                                 breakpoints = self.compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)
1196                              elif self._breakpoint_strategy == 3:
1197                                 breakpoints = self.compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, self._linearize_nonbinary_penalty_terms)                                                                                                   
1198                              else:
1199                                 raise ValueError, "A breakpoint distribution strategy="+str(self._breakpoint_strategy)+" is currently not supported within PH!"
1200
1201                              for i in range(0,len(breakpoints)-1):
1202
1203                                 this_lb = breakpoints[i]
1204                                 this_ub = breakpoints[i+1]
1205
1206                                 piece_constraint_name = "QUAD_PENALTY_PIECE_"+str(i)+"_"+variable_name+str(index)
1207                                 if hasattr(instance, piece_constraint_name) is False:
1208                                    # this is the first time the constraint is being added - add it to the list of PH-specific constraints for this instance.
1209                                    self._instance_augmented_attributes[instance_name].append(piece_constraint_name)                                                                           
1210                                 piece_constraint = Constraint(name=piece_constraint_name)
1211                                 piece_constraint.model = instance
1212                                 piece_expression = self._create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, quad_penalty_term_variable[index])
1213                                 piece_constraint._initialize_constraint(None, piece_expression, instance, piece_constraint_name)
1214                                 setattr(instance, piece_constraint_name, piece_constraint)                                   
1215
1216                           else:
1217
1218                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
1219                   
1220         # strictly speaking, this probably isn't necessary - parameter coefficients won't get
1221         # pre-processed out of the expression tree. however, if the under-the-hood should change,
1222         # we'll be covered.
1223         objective_expression.simplify(instance)
1224         instance._component[Objective][objective_name]._data[None].expr = objective_expression
1225         # if we are linearizing everything, then nothing will appear in the quadratic expression -
1226         # don't add the empty "0.0" expression to the objective. otherwise, the output file won't
1227         # be properly generated.
1228         if quad_expression != 0.0:
1229           instance._component[Objective][objective_name]._quad_subexpr = quad_expression
1230
1231   def iteration_k_solve(self):
1232
1233     if self._verbose is True:
1234        print "------------------------------------------------"       
1235        print "Starting PH iteration " + str(self._current_iteration) + " solves"
1236
1237     # cache the objective values generated by PH for output at the end of this function.
1238     ph_objective_values = {}
1239
1240     solve_start_time = time.time()
1241
1242     # STEP 1: queue up the solves for all scenario sub-problems and
1243     #         grab all of the action handles for the subsequent barrier sync.
1244
1245     action_handles = []
1246     action_handle_instance_map = {}
1247
1248     for scenario in self._scenario_tree._scenarios:     
1249
1250        instance = self._instances[scenario._name]
1251
1252        if self._verbose is True:
1253           print "Queuing solve for scenario=" + scenario._name
1254
1255        # IMPT: You have to re-presolve, as the simple presolver collects the linear terms together. If you
1256        # don't do this, you won't see any chance in the output files as you vary the problem parameters!
1257        # ditto for instance fixing!
1258        instance.presolve()
1259
1260        # once past iteration 0, there is always a feasible solution from which to warm-start.
1261        # however, you might want to disable warm-start when the solver is behaving badly (which does happen).
1262        new_action_handle = None
1263        if self._disable_warmstarts is False:
1264           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, warmstart=True, tee=self._output_solver_log)
1265        else:
1266           new_action_handle = self._solver_manager.queue(instance, opt=self._solver, tee=self._output_solver_log)           
1267
1268        action_handle_instance_map[scenario._name] = new_action_handle
1269
1270        action_handles.append(new_action_handle)
1271
1272     # STEP 2: barrier sync for all scenario sub-problem solves.
1273     if self._verbose is True:           
1274        print "Waiting for scenario sub-problem solves"
1275     self._solver_manager.wait_all(action_handles)
1276     if self._verbose is True:               
1277        print "Scenario sub-problem solves completed"
1278       
1279     solve_end_time = time.time()
1280     self._cumulative_solve_time += (solve_end_time - solve_start_time)
1281
1282     if self._output_times is True:
1283        print "Aggregate sub-problem solve time=%8.2f" % (solve_end_time - solve_start_time)
1284
1285     # STEP 3: Load the results!
1286     for scenario_name, action_handle in action_handle_instance_map.items():
1287
1288        if self._verbose is True:         
1289           print "Successfully processed results for scenario="+scenario_name
1290
1291        instance = self._instances[scenario_name]
1292        results = self._solver_manager.get_results(action_handle)
1293
1294        if len(results.solution) == 0:
1295           raise RuntimeError, "Solve failed for scenario="+scenario_name+"; no solutions generated"
1296
1297        if self._output_solver_results is True:
1298           print "Results for scenario=", scenario_name
1299           results.write(num=1)           
1300
1301        instance.load(results)
1302
1303        if self._verbose is True:                 
1304           print "Successfully loaded solution for scenario="+scenario_name
1305
1306        # we're assuming there is a single solution.
1307        # the "value" attribute is a pre-defined feature of any solution - it is relative to whatever
1308        # objective was selected during optimization, which of course should be the PH objective.
1309        ph_objective_values[instance.name] = float(results.solution(0).objective['f'].value)
1310
1311     if self._verbose is True:
1312        print "Successfully completed PH iteration " + str(self._current_iteration) + " solves - solution statistics:"
1313        print "  Scenario             PH Objective             Cost Objective"
1314        for scenario in self._scenario_tree._scenarios:
1315           instance = self._instances[scenario._name]
1316           for objective_name in instance._component[Objective]:
1317              objective = instance._component[Objective][objective_name]
1318              print "%20s       %18.4f     %14.4f" % (scenario._name, ph_objective_values[scenario._name], 0.0)
1319
1320   def solve(self):
1321
1322      self._solve_start_time = time.time()
1323      self._cumulative_solve_time = 0.0
1324      self._cumulative_xbar_time = 0.0
1325      self._cumulative_weight_time = 0.0
1326
1327      print "Starting PH"
1328
1329      if self._initialized == False:
1330         raise RuntimeError, "PH is not initialized - cannot invoke solve() method"
1331
1332      print "Initiating PH iteration=" + `self._current_iteration`
1333
1334      self.iteration_0_solve()
1335
1336      # update variable statistics prior to any output.
1337      self.update_variable_statistics()     
1338
1339      if self._verbose is True:
1340         print "Variable values following scenario solves:"
1341         self.pprint(False,False,True,False)
1342
1343      # let plugins know if they care.
1344      for plugin in self._ph_plugins:     
1345         plugin.post_iteration_0_solves(self)
1346
1347      # update the fixed variable statistics.
1348      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1349
1350      if self._verbose is True:
1351         print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1352         print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1353
1354      self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1355      print "Convergence metric=%12.4f" % self._converger.lastMetric()
1356
1357      self.update_weights()
1358
1359      # let plugins know if they care.
1360      for plugin in self._ph_plugins:     
1361         plugin.post_iteration_0(self)
1362
1363      # checkpoint if it's time - which it always is after iteration 0,
1364      # if the interval is >= 1!
1365      if (self._checkpoint_interval > 0):
1366         self.checkpoint(0)
1367
1368      # there is an upper bound on the number of iterations to execute -
1369      # the actual bound depends on the converger supplied by the user.
1370      for i in range(1, self._max_iterations+1):
1371
1372         self._current_iteration = self._current_iteration + 1                 
1373
1374         print "Initiating PH iteration=" + `self._current_iteration`         
1375
1376         if self._verbose is True:
1377            print "Variable averages and weights prior to scenario solves:"
1378            self.pprint(True,True,False,False)
1379
1380         # with the introduction of piecewise linearization, the form of the
1381         # penalty-weighted objective is no longer fixed. thus, we need to
1382         # create the objectives each PH iteration.
1383         self.form_iteration_k_objectives()           
1384
1385         self.iteration_k_solve()
1386
1387         # update variable statistics prior to any output.
1388         self.update_variable_statistics()
1389         
1390         if self._verbose is True:
1391            print "Variable values following scenario solves:"
1392            self.pprint(False,False,True,False)
1393
1394         # we don't technically have to do this at the last iteration,
1395         # but with checkpointing and re-starts, you're never sure
1396         # when you're executing the last iteration.
1397         self.update_weights()
1398
1399         # let plugins know if they care.
1400         for plugin in self._ph_plugins:
1401            plugin.post_iteration_k_solves(self)
1402
1403         # update the fixed variable statistics.
1404         (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()               
1405
1406         if self._verbose is True:
1407            print "Number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1408            print "Number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"
1409
1410         # let plugins know if they care.
1411         for plugin in self._ph_plugins:
1412            plugin.post_iteration_k(self)
1413
1414         # at this point, all the real work of an iteration is complete.
1415
1416         # checkpoint if it's time.
1417         if (self._checkpoint_interval > 0) and (i % self._checkpoint_interval is 0):
1418            self.checkpoint(i)
1419
1420         # check for early termination.
1421         self._converger.update(self._current_iteration, self, self._scenario_tree, self._instances)
1422         print "Convergence metric=%12.4f" % self._converger.lastMetric()
1423
1424         if self._converger.isConverged(self) is True:
1425            if self._total_discrete_vars == 0:
1426               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)
1427            else:
1428               print "PH converged - convergence metric is below threshold="+str(self._converger._convergence_threshold)+" or all discrete variables are fixed"               
1429            break
1430
1431         # if we're terminating due to exceeding the maximum iteration count, print a message
1432         # indicating so - otherwise, you get a quiet, information-free output trace.
1433         if i == self._max_iterations:
1434            print "Halting PH - reached maximal iteration count="+str(self._max_iterations)
1435
1436      if self._verbose is True:
1437         print "Number of discrete variables fixed before final plugin calls="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1438         print "Number of continuous variables fixed before final plugin calls="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"           
1439
1440      # let plugins know if they care. do this before
1441      # the final solution / statistics output, as the plugins
1442      # might do some final tweaking.
1443      for plugin in self._ph_plugins:
1444         plugin.post_ph_execution(self)
1445
1446      # update the fixed variable statistics - the plugins might have done something.
1447      (self._total_fixed_discrete_vars,self._total_fixed_continuous_vars) = self.compute_fixed_variable_counts()                       
1448
1449      print "PH complete"
1450
1451      print "Convergence history:"
1452      self._converger.pprint()
1453
1454      print "Final number of discrete variables fixed="+str(self._total_fixed_discrete_vars)+" (total="+str(self._total_discrete_vars)+")"
1455      print "Final number of continuous variables fixed="+str(self._total_fixed_continuous_vars)+" (total="+str(self._total_continuous_vars)+")"         
1456
1457      print "Final variable values:"
1458      self.pprint(False,False,True,True)         
1459
1460      print "Final costs:"
1461      self._scenario_tree.pprintCosts(self._instances)
1462
1463      self._solve_end_time = time.time()
1464
1465      if (self._verbose is True) and (self._output_times is True):
1466         print "Overall run-time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1467
1468      # cleanup the scenario instances for post-processing - ideally, we want to leave them in
1469      # their original state, minus all the PH-specific stuff. we don't do all cleanup (leaving
1470      # things like rhos, etc), but we do clean up constraints, as that really hoses up the ef writer.
1471      self._cleanup_scenario_instances()
1472
1473   #
1474   # prints a summary of all collected time statistics
1475   #
1476   def print_time_stats(self):
1477
1478      print "PH run-time statistics (user):"
1479
1480      print "Initialization time=  %8.2f seconds" % (self._init_end_time - self._init_start_time)
1481      print "Overall solve time=   %8.2f seconds" % (self._solve_end_time - self._solve_start_time)
1482      print "Scenario solve time=  %8.2f seconds" % self._cumulative_solve_time
1483      print "Average update time=  %8.2f seconds" % self._cumulative_xbar_time
1484      print "Weight update time=   %8.2f seconds" % self._cumulative_weight_time
1485
1486   #
1487   # a utility to determine whether to output weight / average / etc. information for
1488   # a variable/node combination. when the printing is moved into a callback/plugin,
1489   # this routine will go there. for now, we don't dive down into the node resolution -
1490   # just the variable/stage.
1491   #
1492   def should_print(self, stage, variable, variable_indices):
1493
1494      if self._output_continuous_variable_stats is False:
1495
1496         variable_type = variable.domain         
1497
1498         if (isinstance(variable_type, IntegerSet) is False) and (isinstance(variable_type, BooleanSet) is False):
1499
1500            return False
1501
1502      return True
1503     
1504   #
1505   # pretty-prints the state of the current variable averages, weights, and values.
1506   # inputs are booleans indicating which components should be output.
1507   #
1508   def pprint(self, output_averages, output_weights, output_values, output_fixed):
1509
1510      if self._initialized is False:
1511         raise RuntimeError, "PH is not initialized - cannot invoke pprint() method"         
1512     
1513      # print tree nodes and associated variable/xbar/ph information in stage-order
1514      # we don't blend in the last stage, so we don't current care about printing the associated information.
1515      for stage in self._scenario_tree._stages[:-1]:
1516
1517         print "\tStage=" + stage._name
1518
1519         num_outputs_this_stage = 0 # tracks the number of outputs on a per-index basis.
1520
1521         for (variable, index_template, variable_indices) in stage._variables:
1522
1523            variable_name = variable.name
1524
1525            if self.should_print(stage, variable, variable_indices) is True:
1526
1527               num_outputs_this_variable = 0 # track, so we don't output the variable names unless there is an entry to report.
1528
1529               for index in variable_indices:               
1530
1531                  weight_parameter_name = "PHWEIGHT_"+variable_name
1532
1533                  num_outputs_this_index = 0 # track, so we don't output the variable index more than once.
1534
1535                  for tree_node in stage._tree_nodes:                 
1536
1537                     # determine if the variable/index pair is used across the set of scenarios (technically,
1538                     # it should be good enough to check one scenario). ditto for "fixed" status. fixed does
1539                     # imply unused (see note below), but we care about the fixed status when outputting
1540                     # final solutions.
1541
1542                     is_used = True # should be consistent across scenarios, so one "unused" flags as invalid.
1543                     is_fixed = False
1544
1545                     for scenario in tree_node._scenarios:
1546                        instance = self._instances[scenario._name]
1547                        variable_value = getattr(instance,variable_name)[index]
1548                        if variable_value.status == VarStatus.unused:
1549                           is_used = False
1550                        if variable_value.fixed is True:
1551                           is_fixed = True
1552 
1553                     # IMPT: this is far from obvious, but variables that are fixed will - because
1554                     #       presolve will identify them as constants and eliminate them from all
1555                     #       expressions - be flagged as "unused" and therefore not output.
1556
1557                     if ((output_fixed is True) and (is_fixed is True)) or (is_used is True):
1558
1559                           minimum_value = tree_node._minimums[variable_name][index]()
1560                           maximum_value = tree_node._maximums[variable_name][index]()
1561
1562                           num_outputs_this_stage = num_outputs_this_stage + 1                           
1563                           num_outputs_this_variable = num_outputs_this_variable + 1
1564                           num_outputs_this_index = num_outputs_this_index + 1
1565
1566                           if num_outputs_this_variable == 1:
1567                              print "\t\tVariable=",variable_name
1568
1569                           if num_outputs_this_index == 1:
1570                              print "\t\t\tIndex:", indexToString(index)                             
1571
1572                           print "\t\t\t\tTree Node=",tree_node._name,"\t\t (Scenarios: ",                             
1573                           for scenario in tree_node._scenarios:
1574                              print scenario._name," ",
1575                              if scenario == tree_node._scenarios[-1]:
1576                                 print ")"
1577                           
1578                           if output_values is True:
1579                              average_value = tree_node._averages[variable_name][index]()
1580                              print "\t\t\t\tValues: ",                       
1581                              for scenario in tree_node._scenarios:
1582                                 instance = self._instances[scenario._name]
1583                                 this_value = getattr(instance,variable_name)[index].value
1584                                 print "%12.4f" % this_value,
1585                                 if scenario == tree_node._scenarios[-1]:
1586                                    print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1587                                    print "    Avg=%12.4f" % (average_value),
1588                                    print ""
1589                           if output_weights:
1590                              print "\t\t\t\tWeights: ",
1591                              for scenario in tree_node._scenarios:
1592                                 instance = self._instances[scenario._name]
1593                                 print "%12.4f" % getattr(instance,weight_parameter_name)[index].value,
1594                                 if scenario == tree_node._scenarios[-1]:
1595                                    print ""
1596
1597                           if output_averages:
1598                              print "\t\t\t\tAverage: %12.4f" % (tree_node._averages[variable_name][index].value)
1599
1600         if num_outputs_this_stage == 0:
1601            print "\t\tNo non-converged variables in stage"
1602
1603         # cost variables aren't blended, so go through the gory computation of min/max/avg.
1604         # we currently always print these.
1605         cost_variable_name = stage._cost_variable[0].name
1606         cost_variable_index = stage._cost_variable[1]
1607         if cost_variable_index is None:
1608            print "\t\tCost Variable=" + cost_variable_name
1609         else:
1610            print "\t\tCost Variable=" + cost_variable_name + indexToString(cost_variable_index)           
1611         for tree_node in stage._tree_nodes:
1612            print "\t\t\tTree Node=" + tree_node._name + "\t\t (Scenarios: ",
1613            for scenario in tree_node._scenarios:
1614               print scenario._name," ",
1615               if scenario == tree_node._scenarios[-1]:
1616                  print ")"
1617            maximum_value = 0.0
1618            minimum_value = 0.0
1619            sum_values = 0.0
1620            num_values = 0
1621            first_time = True
1622            print "\t\t\tValues: ",                       
1623            for scenario in tree_node._scenarios:
1624                instance = self._instances[scenario._name]
1625                this_value = getattr(instance,cost_variable_name)[cost_variable_index].value
1626                print "%12.4f" % this_value,
1627                num_values += 1
1628                sum_values += this_value
1629                if first_time is True:
1630                   first_time = False
1631                   maximum_value = this_value
1632                   minimum_value = this_value
1633                else:
1634                   if this_value > maximum_value:
1635                      maximum_value = this_value
1636                   if this_value < minimum_value:
1637                      minimum_value = this_value
1638                if scenario == tree_node._scenarios[-1]:
1639                   print "    Max-Min=%12.4f" % (maximum_value-minimum_value),
1640                   print "    Avg=%12.4f" % (sum_values/num_values),
1641                   print ""
1642           
Note: See TracBrowser for help on using the repository browser.