source: coopr.pysp/trunk/coopr/pysp/phobjective.py @ 2405

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

Restructuring of PySP to facilitate full implementation of PH solver servers (penalty objectives in particular).

  • Property svn:executable set to *
File size: 19.0 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
4#  Copyright (c) 2010 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
11# this module contains various utilities for creating PH weighted penalty
12# objectives, e.g., through quadratic or linearized penalty terms.
13
14from math import fabs, log, exp
15from coopr.pyomo import *
16
17# IMPT: In general, the breakpoint computation codes can return a 2-list even if the lb equals
18#       the ub. This case happens quite often in real models (although typically lb=xvag=ub).
19#       See the code for constructing the pieces on how this case is handled in the linearization.
20
21#
22# routine to compute linearization breakpoints uniformly between the bounds and the mean.
23#
24
25def compute_uniform_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints_per_side, tolerance):
26   
27   breakpoints = []
28
29   # add the lower bound - the first breakpoint.
30   breakpoints.append(lb)
31
32   # determine the breakpoints to the left of the mean.
33   left_step = (xavg - lb) / num_breakpoints_per_side
34   current_x = lb
35   for i in range(1,num_breakpoints_per_side+1):
36      this_lb = current_x
37      this_ub = current_x+left_step
38      if (fabs(this_lb-lb) > tolerance) and (fabs(this_lb-xavg) > tolerance):
39         breakpoints.append(this_lb)
40      current_x += left_step           
41
42   # add the mean - it's always a breakpoint. unless!
43   # the lb or ub and the avg are the same.
44   if (fabs(lb-xavg) > tolerance) and (fabs(ub-xavg) > tolerance):
45      breakpoints.append(xavg)
46
47   # determine the breakpoints to the right of the mean.
48   right_step = (ub - xavg) / num_breakpoints_per_side
49   current_x = xavg
50   for i in range(1,num_breakpoints_per_side+1):
51      this_lb = current_x
52      this_ub = current_x+right_step
53      if (fabs(this_ub-xavg) > tolerance) and (fabs(this_ub-ub) > tolerance):         
54         breakpoints.append(this_ub)
55      current_x += right_step
56
57   # add the upper bound - the last breakpoint.
58   # the upper bound should always be different than the lower bound (I say with some
59   # hesitation - it really depends on what plugins are doing to modify the bounds dynamically).
60   breakpoints.append(ub)
61
62   return breakpoints
63
64#
65# routine to compute linearization breakpoints uniformly between the current node min/max bounds.
66#   
67
68def compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints, tolerance):
69
70   breakpoints = []
71
72   # add the lower bound - the first breakpoint.
73   breakpoints.append(lb)
74
75   # add the node-min - the second breakpoint. but only if it is different than the lower bound and the mean.
76   if (fabs(node_min-lb) > tolerance) and (fabs(node_min-xavg) > tolerance):     
77      breakpoints.append(node_min)
78
79   step = (node_max - node_min) / num_breakpoints
80   current_x = node_min
81   for i in range(1,num_breakpoints+1):
82      this_lb = current_x
83      this_ub = current_x+step
84      if (fabs(this_lb-node_min) > tolerance) and (fabs(this_lb-node_max) > tolerance) and (fabs(this_lb-xavg) > tolerance):
85         breakpoints.append(this_lb)
86      current_x += step           
87
88   # add the node-max - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
89   if (fabs(node_max-ub) > tolerance) and (fabs(node_max-xavg) > tolerance):           
90      breakpoints.append(node_max)     
91
92   # add the upper bound - the last breakpoint.
93   breakpoints.append(ub)
94
95   # add the mean - it's always a breakpoint. unless! -
96   # it happens to be equal to (within tolerance) the lower or upper bounds.
97   # sort to insert it in the correct place.
98   if (fabs(xavg - lb) > tolerance) and (fabs(xavg - ub) > tolerance):
99      breakpoints.append(xavg)
100   breakpoints.sort()
101
102   return breakpoints
103
104#
105# routine to compute linearization breakpoints using "Woodruff" relaxation of the compute_uniform_between_nodestat_breakpoints.
106#   
107
108def compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints, tolerance):
109
110   breakpoints = []
111
112   # add the lower bound - the first breakpoint.
113   breakpoints.append(lb)
114
115   # be either three "differences" from the mean, or else "halfway to the bound", whichever is closer to the mean.
116   left = max(xavg - 3.0 * (xavg - node_min), xavg - 0.5 * (xavg - lb))
117   right = min(xavg + 3.0 * (node_max - xavg), xavg + 0.5 * (ub - xavg))     
118
119   # add the left bound - the second breakpoint. but only if it is different than the lower bound and the mean.
120   if (fabs(left-lb) > tolerance) and (fabs(left-xavg) > tolerance):     
121      breakpoints.append(left)
122
123   step = (right - left) / num_breakpoints
124   current_x = left
125   for i in range(1,num_breakpoints+1):
126      this_lb = current_x
127      this_ub = current_x+step
128      if (fabs(this_lb-left) > tolerance) and (fabs(this_lb-right) > tolerance) and (fabs(this_lb-xavg) > tolerance):
129         breakpoints.append(this_lb)
130      current_x += step           
131
132   # add the right bound - the second-to-last breakpoint. but only if it is different than the upper bound and the mean.
133   if (fabs(right-ub) > tolerance) and (fabs(right-xavg) > tolerance):           
134      breakpoints.append(right)     
135
136   # add the upper bound - the last breakpoint.
137   breakpoints.append(ub)
138
139   # add the mean - it's always a breakpoint.
140   # sort to insert it in the correct place.
141   breakpoints.append(xavg)
142   breakpoints.sort()
143
144   return breakpoints
145
146#
147# routine to compute linearization breakpoints based on an exponential distribution from the mean in each direction.
148#   
149
150def compute_exponential_from_mean_breakpoints(lb, node_min, xavg, node_max, ub, num_breakpoints_per_side, tolerance):
151
152   breakpoints = []
153
154   # add the lower bound - the first breakpoint.
155   breakpoints.append(lb)
156
157   # determine the breakpoints to the left of the mean.
158   left_delta = xavg - lb
159   base = exp(log(left_delta) / num_breakpoints_per_side)
160   current_offset = base
161   for i in range(1,num_breakpoints_per_side+1):
162      current_x = xavg - current_offset
163      if (fabs(current_x-lb) > tolerance) and (fabs(current_x-xavg) > tolerance):
164         breakpoints.append(current_x)
165      current_offset *= base
166
167   # add the mean - it's always a breakpoint.
168   breakpoints.append(xavg)
169
170   # determine the breakpoints to the right of the mean.
171   right_delta = ub - xavg
172   base = exp(log(right_delta) / num_breakpoints_per_side)
173   current_offset = base
174   for i in range(1,num_breakpoints_per_side+1):
175      current_x = xavg + current_offset
176      if (fabs(current_x-xavg) > tolerance) and (fabs(current_x-ub) > tolerance):         
177         breakpoints.append(current_x)
178      current_offset *= base
179
180   # add the upper bound - the last breakpoint.
181   breakpoints.append(ub)         
182
183   return breakpoints     
184
185#
186# a utility to create piece-wise linear constraint expressions for a given variable, for
187# use in constructing the augmented (penalized) PH objective. lb and ub are the bounds
188# on this piece, variable is the actual instance variable, and average is the instance
189# parameter specifying the average of this variable across instances sharing passing
190# through a common tree node. lb and ub are floats.
191# IMPT: There are cases where lb=ub, in which case the slope is 0 and the intercept
192#       is simply the penalty at the lower(or upper) bound.
193#
194
195def create_piecewise_constraint_expression(lb, ub, instance_variable, variable_average, quad_variable, tolerance):
196
197   penalty_at_lb = (lb - variable_average()) * (lb - variable_average())
198   penalty_at_ub = (ub - variable_average()) * (ub - variable_average())
199   slope = None
200   if fabs(ub-lb) > tolerance:
201      slope = (penalty_at_ub - penalty_at_lb) / (ub - lb)
202   else:
203      slope = 0.0
204   intercept = penalty_at_lb - slope * lb
205   expression = (0.0, quad_variable - slope * instance_variable - intercept, None)
206
207   return expression
208
209#
210# form the baseline objective. really just a wrapper around a clone
211# operation at this point.
212#
213
214def form_standard_objective(instance_name, instance, original_objective_expression, scenario_tree):
215
216   objective_name = instance.active_components(Objective).keys()[0]         
217   objective = instance.active_components(Objective)[objective_name]
218   # clone the objective, because we're going to augment (repeatedly) the original objective.
219   objective_expression = original_objective_expression.clone() 
220   objective_expression.simplify(instance)
221   instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
222   instance.active_components(Objective)[objective_name]._quad_subexpr = None
223
224#
225# form the penalized PH objective, guided by various options.
226# returns the list of components added to the instance, e.g.,
227# in the case of constraints required to implement linearization.
228#
229
230def form_ph_objective(instance_name, instance, original_objective_expression, scenario_tree, \
231                      linearize_nonbinary_penalty_terms, drop_proximal_terms, \
232                      retain_quadratic_binary_terms, breakpoint_strategy, \
233                      tolerance):
234
235   new_instance_attributes = []
236
237   objective_name = instance.active_components(Objective).keys()[0]         
238   objective = instance.active_components(Objective)[objective_name]
239   # clone the objective, because we're going to augment (repeatedly) the original objective.
240   objective_expression = original_objective_expression.clone() 
241   # the quadratic expression is really treated as just a list - eventually should be treated as a full expression.
242   quad_expression = 0.0
243
244   # for each blended variable (i.e., those not appearing in the final stage),
245   # add the linear and quadratic penalty terms to the objective.
246   
247   for stage in scenario_tree._stages[:-1]: # skip the last stage, as no blending occurs
248
249      variable_tree_node = None
250      for node in stage._tree_nodes:
251         for scenario in node._scenarios:
252            if scenario._name == instance_name:
253               variable_tree_node = node
254               break
255
256      for (reference_variable, index_template, variable_indices) in stage._variables:
257
258         variable_name = reference_variable.name
259         variable_type = reference_variable.domain
260
261         w_parameter_name = "PHWEIGHT_"+variable_name
262         w_parameter = instance.active_components(Param)[w_parameter_name]
263           
264         average_parameter_name = "PHAVG_"+variable_name
265         average_parameter = instance.active_components(Param)[average_parameter_name]
266
267         rho_parameter_name = "PHRHO_"+variable_name
268         rho_parameter = instance.active_components(Param)[rho_parameter_name]
269
270         blend_parameter_name = "PHBLEND_"+variable_name
271         blend_parameter = instance.active_components(Param)[blend_parameter_name]
272
273         node_min_parameter = variable_tree_node._minimums[variable_name]
274         node_max_parameter = variable_tree_node._maximums[variable_name]               
275
276         quad_penalty_term_variable = None
277         if linearize_nonbinary_penalty_terms > 0:
278            quad_penalty_term_variable_name = "PHQUADPENALTY_"+variable_name
279            quad_penalty_term_variable = instance.active_components(Var)[quad_penalty_term_variable_name]
280
281         instance_variable = instance.active_components(Var)[variable_name]
282
283         for index in variable_indices:
284
285            if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
286
287               # add the linear (w-weighted) term is a consistent fashion, independent of variable type.
288               # if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes for reasons currently unknown.
289               objective_expression += (w_parameter[index] * instance_variable[index])
290
291               # there are some cases in which a user may want to avoid the proximal term completely.
292               # it's of course only a good idea when there are at least bounds (both lb and ub) on
293               # the variables to be blended.
294               if drop_proximal_terms is False:
295
296                  # deal with binaries
297                  if isinstance(variable_type, BooleanSet) is True:
298
299                     if retain_quadratic_binary_terms is False:
300                        # this rather ugly form of the linearized quadratic expression term is required
301                        # due to a pyomo bug - the form (rho/2) * (x+y+z) chokes in presolve when distributing
302                        # over the sum.
303                        new_term = (blend_parameter[index] * rho_parameter[index] / 2.0 * instance_variable[index]) - \
304                                   (blend_parameter[index] * rho_parameter[index] * average_parameter[index] * instance_variable[index]) + \
305                                   (blend_parameter[index] * rho_parameter[index] / 2.0 * average_parameter[index] * average_parameter[index])                             
306                        if objective.sense is minimize:
307                           objective_expression += new_term
308                        else:
309                           objective_expression -= new_term                                 
310                     else:
311                        quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
312
313                  # deal with everything else
314                  else:
315
316                     if linearize_nonbinary_penalty_terms > 0:
317
318                        # the variables are easy - just a single entry.
319                        if objective.sense is minimize:
320                           objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
321                        else:
322                           objective_expression -= (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
323                           
324                        # the constraints are somewhat nastier.
325
326                        # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE"
327                        xavg = average_parameter[index]
328                        x = instance_variable[index]
329
330                        lb = None
331                        ub = None
332
333                        if x.lb is None:
334                           raise ValueError, "No lower bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
335                        else:
336                           lb = x.lb()
337
338                        if x.ub is None:
339                           raise ValueError, "No upper bound specified for variable="+variable_name+indexToString(index)+"; required when piece-wise approximating quadratic penalty terms"
340                        else:
341                           ub = x.ub()
342
343                        node_min = node_min_parameter[index]()
344                        node_max = node_max_parameter[index]()
345
346                        # compute the breakpoint sequence according to the specified strategy.
347                        breakpoints = []
348                        if breakpoint_strategy == 0:
349                           breakpoints = compute_uniform_breakpoints(lb, node_min, xavg(), node_max, ub, \
350                                                                     linearize_nonbinary_penalty_terms, \
351                                                                     tolerance)
352                        elif breakpoint_strategy == 1:
353                           breakpoints = compute_uniform_between_nodestat_breakpoints(lb, node_min, xavg(), node_max, ub, \
354                                                                                      linearize_nonbinary_penalty_terms, \
355                                                                                      tolerance)
356                        elif breakpoint_strategy == 2:
357                           breakpoints = compute_uniform_between_woodruff_breakpoints(lb, node_min, xavg(), node_max, ub, \
358                                                                                      linearize_nonbinary_penalty_terms, \
359                                                                                      tolerance)
360                        elif breakpoint_strategy == 3:
361                           breakpoints = compute_exponential_from_mean_breakpoints(lb, node_min, xavg(), node_max, ub, \
362                                                                                   linearize_nonbinary_penalty_terms, \
363                                                                                   tolerance)
364                        else:
365                           raise ValueError, "A breakpoint distribution strategy="+str(breakpoint_strategy)+" is currently not supported within PH!"
366
367                        for i in range(0,len(breakpoints)-1):
368
369                           this_lb = breakpoints[i]
370                           this_ub = breakpoints[i+1]
371
372                           piece_constraint_name = "QUAD_PENALTY_PIECE_"+str(i)+"_"+variable_name+str(index)
373                           if hasattr(instance, piece_constraint_name) is False:
374                              # this is the first time the constraint is being added - add it to the list of PH-specific constraints for this instance.
375                              new_instance_attributes.append(piece_constraint_name)
376
377                           piece_constraint = Constraint(name=piece_constraint_name)
378                           piece_constraint.model = instance
379                           piece_expression = create_piecewise_constraint_expression(this_lb, this_ub, x, xavg, \
380                                                                                     quad_penalty_term_variable[index], \
381                                                                                     tolerance)
382                           piece_constraint.add(None, piece_expression)
383                           setattr(instance, piece_constraint_name, piece_constraint)                                   
384
385                     else:
386
387                        quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
388             
389   # strictly speaking, this probably isn't necessary - parameter coefficients won't get
390   # pre-processed out of the expression tree. however, if the under-the-hood should change,
391   # we'll be covered.
392   objective_expression.simplify(instance)
393   instance.active_components(Objective)[objective_name]._data[None].expr = objective_expression
394   # if we are linearizing everything, then nothing will appear in the quadratic expression -
395   # don't add the empty "0.0" expression to the objective. otherwise, the output file won't
396   # be properly generated.
397   if quad_expression != 0.0:
398     instance.active_components(Objective)[objective_name]._quad_subexpr = quad_expression
399
400   return new_instance_attributes
Note: See TracBrowser for help on using the repository browser.