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

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

Fixed the PH solver server to automatically deregister a delegator object if it already exists in the name server.

Fixed various bugs in the PH solver server relating to quadratic expression output and weight/rho updates. It now replicates CPLEX results across the server network.

  • Property svn:executable set to *
File size: 19.1 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   else:
400     instance.active_components(Objective)[objective_name]._quad_subexpr = None
401
402   return new_instance_attributes
Note: See TracBrowser for help on using the repository browser.