Changeset 1793
 Timestamp:
 Nov 3, 2009 1:23:58 AM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/coopr/pysp/ph.py
r1787 r1793 907 907 # TBD  deal with maximization case!!! 908 908 909 # print "ADDING LINEAR PENALTY TERM TO OBJECTIVE FOR VARIABLE=",instance_variable.name,"INDEX=",index910 909 # 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]) 912 911 913 912 # the constraints are somewhat nastier. 914 913 915 914 # TBD  DEFINE CONSTRAINTS ONTHEFLY?? (INDIVIDUALLY NAMED FOR NOW  CREATE INDEX SETS!)  OR A LEAST AN INDEX SET PER "PIECE" 916 rho = rho_parameter[index]917 915 xavg = average_parameter[index] 918 916 x = instance_variable[index] … … 931 929 # unrealistic bias in terms of computed slope in the piecewise linear 932 930 # 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)) 957 956 958 957 if first_slope() >= 0.0: 959 958 raise RuntimeError, "***ERROR: First piece of linear approximation for quadratic penalty term must have a negative slope! Computed value="+str(first_slope()) 960 959 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)) 973 969 974 970 if second_slope() <= 0.0: 975 971 raise RuntimeError, "***ERROR: Second piece of linear approximation for quadratic penalty term must have a positive slope! Computed value="+str(second_slope()) 976 972 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 985 979 piece_1_expression = (0.0, 986 980 (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) 993 987 ), 994 988 None) … … 998 992 setattr(instance, piece_1_constraint_name, piece_1_constraint) 999 993 1000 # piece_2_expression = (quad_penalty_term_variable[index]  second_slope * x  second_intercept) >= 0.0994 # piece_2_expression = (quad_penalty_term_variable[index]  second_slope * x  second_intercept) >= 0.0 1001 995 piece_2_expression = (0.0, 1002 996 (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)) 1009 1003 ), 1010 1004 None)
Note: See TracChangeset
for help on using the changeset viewer.