Changeset 1784


Ignore:
Timestamp:
Nov 2, 2009 1:26:22 AM (11 years ago)
Author:
jwatson
Message:

Further fixes/enhancements to provide an automatic linearization scheme for quadratic penalty terms involving continuous variables in progressive hedging; includes exploratory changes in presolve to accomodate some more-complex-than-usual expressions.

Location:
trunk/coopr
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/coopr/pyomo/base/constraint.py

    r1768 r1784  
    314314
    315315    def _initialize_constraint(self, val, expr, model, _name):
     316
    316317        #
    317318        # Ignore an 'empty' constraint
  • trunk/coopr/pyomo/base/expr.py

    r1780 r1784  
    411411        tmp._const = self._const
    412412        return tmp
     413    #
     414    # this method contrast with the fixed_value() method.
     415    # the fixed_value() method returns true iff the value is
     416    # an atomic constant.
     417    # this method returns true iff all composite arguments
     418    # in this sum expression are constant, i.e., numeric
     419    # constants or parametrs. the parameter values can of
     420    # course change over time, but at any point in time,
     421    # they are constant. hence, the name.
     422    #
     423    def evaluates_to_constant(self):
     424        for arg in self._args:
     425            if arg.fixed_value() is False:
     426                return False
     427        return True
    413428
    414429    def scale(self, val):
  • trunk/coopr/pyomo/presolve/collect_linear_terms.py

    r1780 r1784  
    6464           #
    6565           elif isinstance(exp,expr._ProductExpression):
    66               print "PROCESSING PRODUCT EXPRESSION:"
    67               exp.pprint()
    6866              c = scale * exp.coef
    6967              ve = v = "0"
    7068              for e in exp._numerator:
    71                  print "PROCESSING NUMERATOR - EXPRESSION:"
    72                  e.pprint()
    73                  print ""                 
    7469                 if e.fixed_value():
    7570                     # at this point, we have checked that the expression has a fixed value.
     
    9590                     ve = e
    9691              for e in exp._denominator:
    97                    print "PROCESSING DENOMINATOR - EXPRESSION:"
    98                    e.pprint()
    99                    print ""
    10092                   if e.fixed_value():
    10193                       c /= e.value
     94                   elif (isinstance(e, expr._SumExpression) is True) and (e.evaluates_to_constant() is True):
     95                      # the following logic should be extended to product expressions and all other types
     96                      # of core expressions; this is experimental, based on issues with logically valid
     97                      # expressions being flagged as "nonlinear, with variables". in particular, the
     98                      # evaluates_to_constant() should be a method of the base class.
     99                      c /= e()
    102100                   else:
    103                        # IF IT'S NOT FIXED, SEE IF THE () VALUE WORKS (is-constant?)
    104                        print "***ERROR***"
    105                        e.pprint()
     101                       print "***Error identified with expression=",e
    106102                       raise ValueError, "Nonlinear expression: found variable in denominator"
    107103              if v in x:
  • trunk/coopr/pysp/ph.py

    r1780 r1784  
    866866                     if (instance_variable[index].status is not VarStatus.unused) and (instance_variable[index].fixed is False):
    867867
     868                        # add the linear (w-weighted) term is a consistent fashion, independent of variable type.
    868869                        # TBD - if maximizing, here is where you would want "-=" - however, if you do this, the collect/simplify process chokes.
    869870                        objective_expression += (w_parameter[index] * instance_variable[index])
     871
     872                        # deal with binaries
    870873
    871874                        if isinstance(variable_type, BooleanSet) is True:
     
    885888                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)
    886889
     890                        # TBD - WARNING - WHAT ABOUT GENERAL INTEGERS? WE CANT LINEARIZE THESE, BUT HOW DO WE DETECT THE SETS?
     891
     892                        # deal with continuous variables.
    887893                        else:
    888894
    889895                           if self._linearize_continuous_variable_penalty_terms is True:
     896
    890897                              # TBD - deal with maximization case!!!
    891898
    892899                              print "ADDING LINEAR PENALTY TERM TO OBJECTIVE FOR VARIABLE=",instance_variable.name,"INDEX=",index                             
    893 
    894900                              # the variables are easy - just a single entry.
    895901                              objective_expression += quad_penalty_term_variable[index]
     
    901907                              xavg = average_parameter[index]
    902908                              x = instance_variable[index]
    903                               # TBD - EXTRACT THESE FROM SOMEWHERE REAL - EITHER THE VARIABLE ITSELF OR THE CROSS-SCENARIO MIN/MAX (PROPERLY ADJUSTED)
    904                               lb = 0.0
    905                               ub = 500.0
     909
     910                              # TBD - MOVE THIS OUTSIDE THE INDEX LOOP - REDUNDANT HERE
     911                              variable_tree_node = None
     912                              for node in stage._tree_nodes:
     913                                 for scenario in node._scenarios:
     914                                    if scenario._name == instance_name:
     915                                       variable_tree_node = node
     916                                       break
     917
     918                              # use the variable min/max as a proxy for the lower and upper bounds.
     919                              # for one, the lower and upper bounds may not be provided. further,
     920                              # given a limited number of segements, they may provide completely
     921                              # unrealistic bias in terms of computed slope in the piecewise linear
     922                              # approximation.
     923                              lb=variable_tree_node._minimums[variable_name][index]
     924                              ub=variable_tree_node._maximums[variable_name][index]
     925
     926                              print "LB=",lb()
    906927                              print "AVERAGE=",xavg()
     928                              print "UB=",ub()
     929
    907930#                              penalty_at_lb = rho / 2.0 * (lb - xavg) * (lb - xavg)
    908931                              penalty_at_lb = (rho / 2.0 * lb * lb)  - \
     
    917940                              print "PENALTY AT LB=", penalty_at_lb()                             
    918941                              print "PENALTY AT UB=", penalty_at_ub()
    919                              
     942
    920943#                              first_slope = -penalty_at_lb / (xavg - lb)
    921                               first_slope = (rho / 2.0 * lb * lb / (xavg - lb)) + \
     944                              first_slope = -(rho / 2.0 * lb * lb / (xavg - lb)) + \
    922945                                            (rho * lb * xavg / (xavg - lb)) - \
    923946                                            (rho / 2.0 * xavg * xavg / (xavg - lb))
    924947
    925 #                              first_intercept = xavg / (xavg - lb) penalty_at_lb
     948                              if first_slope() >= 0.0:
     949                                 raise RuntimeError, "***ERROR: First piece of linear approximation for quadratic penalty term must have a negative slope! Computed value="+str(first_slope())
     950
     951#                              first_intercept = xavg / (xavg - lb) * penalty_at_lb
    926952                              first_intercept = (xavg / (xavg - lb) * rho / 2.0 * lb * lb)  - \
    927953                                                (xavg / (xavg - lb) * rho * lb * xavg) + \
     
    935961                                              (rho * ub * xavg / (ub - xavg)) + \
    936962                                              (rho / 2.0 * xavg * xavg / (ub - xavg))
     963
     964                              if second_slope() <= 0.0:
     965                                 raise RuntimeError, "***ERROR: Second piece of linear approximation for quadratic penalty term must have a positive slope! Computed value="+str(second_slope())
    937966                             
    938967#                              second_intercept = -xavg / (ub - xavg) * penalty_at_ub
     
    954983                              piece_1_constraint_name = "QUAD_PENALTY_PIECE_1_"+variable_name+str(index)
    955984                              piece_1_constraint = Constraint(name=piece_1_constraint_name)
    956                               piece_1_constraint._initialize_constraint(None, piece_1_expression, instance, "")
     985                              piece_1_constraint._initialize_constraint(None, piece_1_expression, instance, piece_1_constraint_name)
    957986                              setattr(instance, piece_1_constraint_name, piece_1_constraint)
    958987
    959988#                              piece_2_expression = (quad_penalty_term_variable[index] - second_slope * x - second_intercept) >= 0.0
    960 # TBD - THE FIRST EXREPSSION IS YIELDING A "FOUND VARIABLE IN DENOMINATOR", WHICH IS FALSE
    961989                              piece_2_expression = (quad_penalty_term_variable[index] - \
    962990                                                    (rho / 2.0 * ub * ub / (ub - xavg) * x) + \
     
    969997                              piece_2_constraint_name = "QUAD_PENALTY_PIECE_2_"+variable_name+str(index)
    970998                              piece_2_constraint = Constraint(name=piece_2_constraint_name)
    971                               piece_2_constraint._initialize_constraint(None, piece_2_expression, instance, "")
     999                              piece_2_constraint._initialize_constraint(None, piece_2_expression, instance, piece_2_constraint_name)
    9721000                              setattr(instance, piece_2_constraint_name, piece_2_constraint)                             
    9731001
    9741002                           else:
     1003
    9751004                              quad_expression += (blend_parameter[index] * rho_parameter[index] * (instance_variable[index] - average_parameter[index]) ** 2)                           
    9761005                   
Note: See TracChangeset for help on using the changeset viewer.