Changeset 1793


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

More improvements/fixes to PH proximal term linearization.

File:
1 edited

Legend:

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

    r1787 r1793  
    907907                                 # TBD - deal with maximization case!!!
    908908
    909    #                              print "ADDING LINEAR PENALTY TERM TO OBJECTIVE FOR VARIABLE=",instance_variable.name,"INDEX=",index                             
    910909                                 # the variables are easy - just a single entry.
    911                                  objective_expression += quad_penalty_term_variable[index]
     910                                 objective_expression += (rho_parameter[index] / 2.0 * quad_penalty_term_variable[index])
    912911
    913912                                 # the constraints are somewhat nastier.
    914913
    915914                                 # TBD - DEFINE CONSTRAINTS ON-THE-FLY?? (INDIVIDUALLY NAMED FOR NOW - CREATE INDEX SETS!) - OR A LEAST AN INDEX SET PER "PIECE"
    916                                  rho = rho_parameter[index]
    917915                                 xavg = average_parameter[index]
    918916                                 x = instance_variable[index]
     
    931929                                 # unrealistic bias in terms of computed slope in the piecewise linear
    932930                                 # approximation.
    933                                  lb=variable_tree_node._minimums[variable_name][index]
    934                                  ub=variable_tree_node._maximums[variable_name][index]
    935 
    936 #                              print "LB=",lb()
    937 #                              print "AVERAGE=",xavg()
    938 #                              print "UB=",ub()
    939 
    940 #                              penalty_at_lb = rho / 2.0 * (lb - xavg) * (lb - xavg)
    941                                  penalty_at_lb = (rho / 2.0 * lb * lb)  - \
    942                                                  (rho * lb * xavg) + \
    943                                                  (rho / 2.0 * xavg * xavg)
    944 
    945 #                              penalty_at_ub = rho / 2.0 * (ub - xavg) * (ub - xavg)
    946                                  penalty_at_ub = (rho / 2.0 * ub * ub) - \
    947                                                  (rho * ub * xavg) + \
    948                                                  (rho / 2.0 * xavg * xavg)
    949 
    950 #                              print "PENALTY AT LB=", penalty_at_lb()                             
    951 #                              print "PENALTY AT UB=", penalty_at_ub()
    952 
    953 #                              first_slope = -penalty_at_lb / (xavg - lb)
    954                                  first_slope = -(rho / 2.0 * lb * lb / (xavg - lb)) + \
    955                                                (rho * lb * xavg / (xavg - lb)) - \
    956                                                (rho / 2.0 * xavg * xavg / (xavg - lb))
     931                                 lb = None
     932                                 ub = None
     933
     934                                 if x.lb is not None:
     935                                    lb = x.lb
     936                                 else:
     937                                    lb=variable_tree_node._minimums[variable_name][index]
     938
     939                                 if x.ub is not None:
     940                                    ub = x.ub
     941                                 else:
     942                                    # TBD - we discussed adding the difference between the ub and lb
     943                                    #       to the ub, in order to give some estimate for the variability.
     944                                    ub=variable_tree_node._maximums[variable_name][index]
     945
     946#                                 penalty_at_lb = (lb - xavg) * (lb - xavg)
     947                                 penalty_at_lb = (lb * lb) - (2.0 * lb * xavg) + (xavg * xavg)
     948
     949#                                 penalty_at_ub = (ub - xavg) * (ub - xavg)
     950                                 penalty_at_ub = (ub * ub) - (2.0 * ub * xavg) + (xavg * xavg)
     951
     952#                                 first_slope = -penalty_at_lb / (xavg - lb)
     953                                 first_slope = -(lb * lb / (xavg - lb)) + \
     954                                               (2.0 * lb * xavg / (xavg - lb)) - \
     955                                               (xavg * xavg / (xavg - lb))
    957956
    958957                                 if first_slope() >= 0.0:
    959958                                    raise RuntimeError, "***ERROR: First piece of linear approximation for quadratic penalty term must have a negative slope! Computed value="+str(first_slope())
    960959
    961 #                              first_intercept = xavg / (xavg - lb) * penalty_at_lb
    962                                  first_intercept = (xavg / (xavg - lb) * rho / 2.0 * lb * lb)  - \
    963                                                    (xavg / (xavg - lb) * rho * lb * xavg) + \
    964                                                    (xavg / (xavg - lb) * rho / 2.0 * xavg * xavg)
    965 
    966 #                              print "FIRST SLOPE=",first_slope()                             
    967 #                              print "FIRST INTERCEPT=",first_intercept()
    968                              
    969 #                              second_slope = penalty_at_ub / (ub - xavg)
    970                                  second_slope =  (rho / 2.0 * ub * ub / (ub - xavg)) - \
    971                                                  (rho * ub * xavg / (ub - xavg)) + \
    972                                                  (rho / 2.0 * xavg * xavg / (ub - xavg))
     960#                                 first_intercept = xavg / (xavg - lb) * penalty_at_lb
     961                                 first_intercept = (xavg / (xavg - lb) * lb * lb)  - \
     962                                                   (xavg / (xavg - lb) * 2.0 * lb * xavg) + \
     963                                                   (xavg / (xavg - lb) * xavg * xavg)
     964
     965#                                 second_slope = penalty_at_ub / (ub - xavg)
     966                                 second_slope =  (ub * ub / (ub - xavg)) - \
     967                                                 (2.0 * ub * xavg / (ub - xavg)) + \
     968                                                 (xavg * xavg / (ub - xavg))
    973969
    974970                                 if second_slope() <= 0.0:
    975971                                    raise RuntimeError, "***ERROR: Second piece of linear approximation for quadratic penalty term must have a positive slope! Computed value="+str(second_slope())
    976972                             
    977 #                              second_intercept = -xavg / (ub - xavg) * penalty_at_ub
    978                                  second_intercept =  (rho / 2.0 * ub * ub * -xavg / (ub - xavg)) - \
    979                                                      (rho * ub * xavg * -xavg / (ub - xavg)) + \
    980                                                      (rho / 2.0 * xavg * xavg * -xavg / (ub - xavg))
    981 #                              print "SECOND SLOPE=",second_slope()
    982 #                              print "SECOND INTERCEPT=",second_intercept()
    983 
    984 #                              piece_1_expression = (quad_penalty_term_variable[index] - first_slope * x - first_intercept) >= 0.0
     973#                                 second_intercept = -xavg / (ub - xavg) * penalty_at_ub
     974                                 second_intercept =  (ub * ub * -xavg / (ub - xavg)) - \
     975                                                     (2.0 * ub * xavg * -xavg / (ub - xavg)) + \
     976                                                     (xavg * xavg * -xavg / (ub - xavg))
     977
     978#                                 piece_1_expression = (quad_penalty_term_variable[index] - first_slope * x - first_intercept) >= 0.0
    985979                                 piece_1_expression =  (0.0,
    986980                                                       (quad_penalty_term_variable[index] - \
    987                                                         -(rho / 2.0 * lb * lb / (xavg - lb) * x) - \
    988                                                         (rho * lb * xavg / (xavg - lb) * x) + \
    989                                                         (rho / 2.0 * xavg * xavg / (xavg - lb) * x) - \
    990                                                         (xavg / (xavg - lb) * rho / 2.0 * lb * lb)  + \
    991                                                         (xavg / (xavg - lb) * rho * lb * xavg) - \
    992                                                         (xavg / (xavg - lb) * rho / 2.0 * xavg * xavg)
     981                                                        -(lb * lb / (xavg - lb) * x) - \
     982                                                        (2.0 * lb * xavg / (xavg - lb) * x) + \
     983                                                        (xavg * xavg / (xavg - lb) * x) - \
     984                                                        (xavg / (xavg - lb) * lb * lb)  + \
     985                                                        (xavg / (xavg - lb) * 2.0 * lb * xavg) - \
     986                                                        (xavg / (xavg - lb) * xavg * xavg)
    993987                                                       ),
    994988                                                       None)
     
    998992                                 setattr(instance, piece_1_constraint_name, piece_1_constraint)
    999993
    1000 #                              piece_2_expression = (quad_penalty_term_variable[index] - second_slope * x - second_intercept) >= 0.0
     994#                                 piece_2_expression = (quad_penalty_term_variable[index] - second_slope * x - second_intercept) >= 0.0
    1001995                                 piece_2_expression = (0.0,
    1002996                                                       (quad_penalty_term_variable[index] - \
    1003                                                        (rho / 2.0 * ub * ub / (ub - xavg) * x) + \
    1004                                                        (rho * ub * xavg / (ub - xavg) * x) - \
    1005                                                        (rho / 2.0 * xavg * xavg / (ub - xavg) *x) - \
    1006                                                        (rho / 2.0 * ub * ub * -xavg / (ub - xavg)) + \
    1007                                                        (rho * ub * xavg * -xavg / (ub - xavg)) - \
    1008                                                        (rho / 2.0 * xavg * xavg * -xavg / (ub - xavg))
     997                                                       (ub * ub / (ub - xavg) * x) + \
     998                                                       (2.0 * ub * xavg / (ub - xavg) * x) - \
     999                                                       (xavg * xavg / (ub - xavg) *x) - \
     1000                                                       (ub * ub * -xavg / (ub - xavg)) + \
     1001                                                       (2.0 * ub * xavg * -xavg / (ub - xavg)) - \
     1002                                                       (xavg * xavg * -xavg / (ub - xavg))
    10091003                                                       ),
    10101004                                                       None)
Note: See TracChangeset for help on using the changeset viewer.