Changeset 1793

Ignore:
Timestamp:
Nov 3, 2009 1:23:58 AM (11 years ago)
Message:

More improvements/fixes to PH proximal term linearization.

File:
1 edited

Unmodified
Added
Removed
• trunk/coopr/pysp/ph.py

 r1787 # TBD - deal with maximization case!!! #                              print "ADDING LINEAR PENALTY TERM TO OBJECTIVE FOR VARIABLE=",instance_variable.name,"INDEX=",index # the variables are easy - just a single entry. objective_expression += quad_penalty_term_variable[index] objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index]) # the constraints are somewhat nastier. # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE" rho = rho_parameter[index] xavg = average_parameter[index] x = instance_variable[index] # unrealistic bias in terms of computed slope in the piecewise linear # approximation. lb=variable_tree_node._minimums[variable_name][index] ub=variable_tree_node._maximums[variable_name][index] #                              print "LB=",lb() #                              print "AVERAGE=",xavg() #                              print "UB=",ub() #                              penalty_at_lb = rho / 2.0 * (lb - xavg) * (lb - xavg) penalty_at_lb = (rho / 2.0 * lb * lb)  - \ (rho * lb * xavg) + \ (rho / 2.0 * xavg * xavg) #                              penalty_at_ub = rho / 2.0 * (ub - xavg) * (ub - xavg) penalty_at_ub = (rho / 2.0 * ub * ub) - \ (rho * ub * xavg) + \ (rho / 2.0 * xavg * xavg) #                              print "PENALTY AT LB=", penalty_at_lb() #                              print "PENALTY AT UB=", penalty_at_ub() #                              first_slope = -penalty_at_lb / (xavg - lb) first_slope = -(rho / 2.0 * lb * lb / (xavg - lb)) + \ (rho * lb * xavg / (xavg - lb)) - \ (rho / 2.0 * xavg * xavg / (xavg - lb)) lb = None ub = None if x.lb is not None: lb = x.lb else: lb=variable_tree_node._minimums[variable_name][index] if x.ub is not None: ub = x.ub else: # TBD - we discussed adding the difference between the ub and lb #       to the ub, in order to give some estimate for the variability. ub=variable_tree_node._maximums[variable_name][index] #                                 penalty_at_lb = (lb - xavg) * (lb - xavg) penalty_at_lb = (lb * lb) - (2.0 * lb * xavg) + (xavg * xavg) #                                 penalty_at_ub = (ub - xavg) * (ub - xavg) penalty_at_ub = (ub * ub) - (2.0 * ub * xavg) + (xavg * xavg) #                                 first_slope = -penalty_at_lb / (xavg - lb) first_slope = -(lb * lb / (xavg - lb)) + \ (2.0 * lb * xavg / (xavg - lb)) - \ (xavg * xavg / (xavg - lb)) if first_slope() >= 0.0: raise RuntimeError, "***ERROR: First piece of linear approximation for quadratic penalty term must have a negative slope! Computed value="+str(first_slope()) #                              first_intercept = xavg / (xavg - lb) * penalty_at_lb first_intercept = (xavg / (xavg - lb) * rho / 2.0 * lb * lb)  - \ (xavg / (xavg - lb) * rho * lb * xavg) + \ (xavg / (xavg - lb) * rho / 2.0 * xavg * xavg) #                              print "FIRST SLOPE=",first_slope() #                              print "FIRST INTERCEPT=",first_intercept() #                              second_slope = penalty_at_ub / (ub - xavg) second_slope =  (rho / 2.0 * ub * ub / (ub - xavg)) - \ (rho * ub * xavg / (ub - xavg)) + \ (rho / 2.0 * xavg * xavg / (ub - xavg)) #                                 first_intercept = xavg / (xavg - lb) * penalty_at_lb first_intercept = (xavg / (xavg - lb) * lb * lb)  - \ (xavg / (xavg - lb) * 2.0 * lb * xavg) + \ (xavg / (xavg - lb) * xavg * xavg) #                                 second_slope = penalty_at_ub / (ub - xavg) second_slope =  (ub * ub / (ub - xavg)) - \ (2.0 * ub * xavg / (ub - xavg)) + \ (xavg * xavg / (ub - xavg)) if second_slope() <= 0.0: raise RuntimeError, "***ERROR: Second piece of linear approximation for quadratic penalty term must have a positive slope! Computed value="+str(second_slope()) #                              second_intercept = -xavg / (ub - xavg) * penalty_at_ub second_intercept =  (rho / 2.0 * ub * ub * -xavg / (ub - xavg)) - \ (rho * ub * xavg * -xavg / (ub - xavg)) + \ (rho / 2.0 * xavg * xavg * -xavg / (ub - xavg)) #                              print "SECOND SLOPE=",second_slope() #                              print "SECOND INTERCEPT=",second_intercept() #                              piece_1_expression = (quad_penalty_term_variable[index] - first_slope * x - first_intercept) >= 0.0 #                                 second_intercept = -xavg / (ub - xavg) * penalty_at_ub second_intercept =  (ub * ub * -xavg / (ub - xavg)) - \ (2.0 * ub * xavg * -xavg / (ub - xavg)) + \ (xavg * xavg * -xavg / (ub - xavg)) #                                 piece_1_expression = (quad_penalty_term_variable[index] - first_slope * x - first_intercept) >= 0.0 piece_1_expression =  (0.0, (quad_penalty_term_variable[index] - \ -(rho / 2.0 * lb * lb / (xavg - lb) * x) - \ (rho * lb * xavg / (xavg - lb) * x) + \ (rho / 2.0 * xavg * xavg / (xavg - lb) * x) - \ (xavg / (xavg - lb) * rho / 2.0 * lb * lb)  + \ (xavg / (xavg - lb) * rho * lb * xavg) - \ (xavg / (xavg - lb) * rho / 2.0 * xavg * xavg) -(lb * lb / (xavg - lb) * x) - \ (2.0 * lb * xavg / (xavg - lb) * x) + \ (xavg * xavg / (xavg - lb) * x) - \ (xavg / (xavg - lb) * lb * lb)  + \ (xavg / (xavg - lb) * 2.0 * lb * xavg) - \ (xavg / (xavg - lb) * xavg * xavg) ), None) setattr(instance, piece_1_constraint_name, piece_1_constraint) #                              piece_2_expression = (quad_penalty_term_variable[index] - second_slope * x - second_intercept) >= 0.0 #                                 piece_2_expression = (quad_penalty_term_variable[index] - second_slope * x - second_intercept) >= 0.0 piece_2_expression = (0.0, (quad_penalty_term_variable[index] - \ (rho / 2.0 * ub * ub / (ub - xavg) * x) + \ (rho * ub * xavg / (ub - xavg) * x) - \ (rho / 2.0 * xavg * xavg / (ub - xavg) *x) - \ (rho / 2.0 * ub * ub * -xavg / (ub - xavg)) + \ (rho * ub * xavg * -xavg / (ub - xavg)) - \ (rho / 2.0 * xavg * xavg * -xavg / (ub - xavg)) (ub * ub / (ub - xavg) * x) + \ (2.0 * ub * xavg / (ub - xavg) * x) - \ (xavg * xavg / (ub - xavg) *x) - \ (ub * ub * -xavg / (ub - xavg)) + \ (2.0 * ub * xavg * -xavg / (ub - xavg)) - \ (xavg * xavg * -xavg / (ub - xavg)) ), None)
Note: See TracChangeset for help on using the changeset viewer.