Changeset 2066


Ignore:
Timestamp:
Dec 29, 2009 10:54:59 PM (10 years ago)
Author:
jwatson
Message:

Removed some TBD notes in the Benders example, and checked in the full benders script (a debug version was committed initially, that only dealt with iteration 1). This version replicates, to within precision, the AMPL execution. Should terminate in 5 master iterations.

Location:
coopr.pyomo/trunk/examples/pyomo/benders
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • coopr.pyomo/trunk/examples/pyomo/benders/base_subproblem.dat

    r1996 r2066  
    1010param prodcost :=  bands 10    coils  11 ;
    1111param invcost  :=  bands  2.5  coils   3 ;
     12
     13param inv1 := bands 0 coils 0 ;
    1214
    1315param market:    1     2     3     4 :=
  • coopr.pyomo/trunk/examples/pyomo/benders/high_subproblem.dat

    r1996 r2066  
    1010param prodcost :=  bands 10    coils  11 ;
    1111param invcost  :=  bands  2.5  coils   3 ;
     12
     13param inv1 := bands 0 coils 0 ;
    1214
    1315param market:    1     2     3     4 :=
  • coopr.pyomo/trunk/examples/pyomo/benders/low_subproblem.dat

    r1996 r2066  
    1010param prodcost :=  bands 10    coils  11 ;
    1111param invcost  :=  bands  2.5  coils   3 ;
     12
     13param inv1 := bands 0 coils 0 ;
    1214
    1315param market:    1     2     3     4 :=
  • coopr.pyomo/trunk/examples/pyomo/benders/master.py

    r2032 r2066  
    5555model.prob = Param(model.SCEN, validate=prob_validator)
    5656
    57 # TBD - we could/should add in a check for the sum of probabilities equal to 1.
    58 
    5957##############################################################################
    6058
     
    9795model.Cut_Defn = Constraint(model.CUTS)
    9896
    99 #maximize Expected_Profit:
    100 #   sum {s in SCEN} prob[s] * 
    101 #     sum {p in PROD} (revenue[p,1,s]*Sell1[p] -
    102 #        prodcost[p]*Make1[p] - invcost[p]*Inv1[p]) +
    103 #   Min_Stage2_Profit;
    104 
    10597def expected_profit_rule(model):
    106    # TBD - still need scenario probability in here!!!
    10798   return sum([model.prob[s] * model.revenue[p, 1, s] * model.Sell1[p] - \
    10899               model.prob[s] * model.prodcost[p] * model.Make1[p] - \
     
    112103
    113104model.Expected_Profit = Objective(rule=expected_profit_rule, sense=maximize)
    114 
    115 
  • coopr.pyomo/trunk/examples/pyomo/benders/runbenders

    r2051 r2066  
    4040GAP = float("Inf")
    4141
    42 max_iterations = 1 # TBD - 50 is what is in the AMPL script
     42max_iterations = 50
    4343
    4444# the main benders loop.
    4545for i in range(1,max_iterations+1):
    4646
    47    print "Iteration=",i
     47   print "\nIteration=",i
    4848
    4949   # solve the subproblems.
    50    print "SOLVING SUB-PROBLEMS"
    5150   solve_all_instances(solver_manager, solver, sub_insts)
    5251   for instance in sub_insts:
    53       print "PROFIT FOR SCENARIO="+instance.name+"="+str(instance.Exp_Stage2_Profit())     
     52      print "Profit for scenario="+instance.name+" is "+str(instance.Exp_Stage2_Profit())     
    5453   print ""
    5554         
     
    5857   for s in range(1, len(sub_insts)+1):
    5958      inst = sub_insts[s-1]
    60       print "Adding cuts for sub-problem=",inst.name
    6159      for t in mstr_inst.TWOPLUSWEEKS:
    6260         mstr_inst.time_price[t, s, i] = inst.Time[t].dual
    63          print "T=",t," Time[t].dual=",inst.Time[t].dual
    6461      for p in mstr_inst.PROD:
    6562         mstr_inst.bal2_price[p, s, i] = inst.Balance2[p].dual
    66          print "P=",p," Balance2[p].dual=",inst.Balance2[p].dual
    6763      for p in mstr_inst.PROD:
    6864         for t in mstr_inst.TWOPLUSWEEKS:
    6965            mstr_inst.sell_lim_price[p, t, s, i] = inst.Sell[p, t].urc
    70             print "Sell["+str(p)+","+str(t)+"].urc="+str(inst.Sell[p,t].urc)
    71       print ""
    7266
    7367   # add the master cut.
     
    8478
    8579   newGAP = mstr_inst.Min_Stage2_Profit.value - Exp_Stage2_Profit
    86    print "NEWGAP = ",newGAP
     80   print "New gap=",newGAP,"\n"
    8781
    8882   if newGAP > 0.00001:
    8983      GAP = min(GAP, newGAP)
    90       print "GAP=",GAP
    9184   else:
    9285      break
    9386
    9487   # re-solve the master and update the subproblem inv1 values.
    95    print "RE-SOLVING THE MASTER PROBLEM!"
    9688   mstr_inst.presolve() # presolve to account for the added cuts!
    97 #   mstr_inst.pprint()
    9889   solve_all_instances(solver_manager, solver, [mstr_inst])
    9990
    100 print "ALL DONE!"
     91   print "Master expected profit="+str(mstr_inst.Expected_Profit())
     92
     93   for instance in sub_insts:
     94      for p in mstr_inst.PROD:
     95         # the master inventory values might be slightly
     96         # less than 0 (within tolerance); threshold here.
     97         instance.inv1[p] = max(mstr_inst.Inv1[p](),0.0)
     98      instance.presolve() # re-generate constraints based on new parameter data
     99
     100print "Benders converged!"
     101
     102print "\nConverged master solution values:"
     103for p in mstr_inst.PROD:
     104   print "Make1["+p+"]="+str(mstr_inst.Make1[p]())
     105   print "Sell1["+p+"]="+str(mstr_inst.Sell1[p]())   
     106   print "Inv1["+p+"]="+str(mstr_inst.Inv1[p]())
  • coopr.pyomo/trunk/examples/pyomo/benders/subproblem.py

    r2023 r2066  
    4949model.prob = Param(validate=unit_interval_validate)
    5050
    51 # inventory at end of first period
    52 # NOTE: this parameter will generally be set by the master problem.
    53 #       however, we need to provide a default value because
    54 #       it is referenced in constraint Balance2. If you
    55 #       don't do this, the index set won't be instantiated,
    56 #       and the model instantiation process will fail.
    57 model.inv1 = Param(model.PROD, within=NonNegativeReals, default=0.0)
     51# inventory at end of first period.
     52model.inv1 = Param(model.PROD, within=NonNegativeReals)
    5853
    5954# tons produced
Note: See TracChangeset for help on using the changeset viewer.