source: coopr.pyomo/trunk/examples/pyomo/benders/runbenders @ 5776

Last change on this file since 5776 was 5776, checked in by jwatson, 8 years ago

Upgrading Benders example to be compliant with recent "attr"->"suffix" change in naming convention. Also eliminated the price validator in master.py, which was bogus to begin with - and was not being properly invoked.

  • Property svn:executable set to *
File size: 3.9 KB
Line 
1#! /usr/bin/env python
2
3#
4# this python script is a pyomo-centric translation of the AMPL
5# script found at: http://www.ampl.com/NEW/LOOP2/stoch2.run
6#
7
8# Python imports
9import sys
10from pyutilib.misc import import_file
11from coopr.opt.base import SolverFactory
12from coopr.opt.parallel import SolverManagerFactory
13from coopr.opt.parallel.manager import solve_all_instances
14from coopr.pyomo import *
15
16# initialize the master instance.
17mstr_mdl = import_file("master.py").model
18mstr_inst = mstr_mdl.create("master.dat")
19
20# initialize the sub-problem instances.
21sb_mdl = import_file("subproblem.py").model
22sub_insts = [] # a Python list
23sub_insts.append(sb_mdl.create(name="Base Sub-Problem", \
24                               filename="base_subproblem.dat"))
25sub_insts.append(sb_mdl.create(name="Low Sub-Problem", \
26                               filename="low_subproblem.dat"))
27sub_insts.append(sb_mdl.create(name="High Sub-Problem", \
28                               filename="high_subproblem.dat"))
29
30# initialize the solver / solver manager.
31solver = SolverFactory("cplex")
32if solver is None:
33    print "A CPLEX solver is not available on this machine."
34    sys.exit(1)
35solver_manager = SolverManagerFactory("serial")
36
37# miscellaneous initialization.
38mstr_inst.Min_Stage2_Profit = float("Inf")
39
40GAP = float("Inf")
41
42max_iterations = 50 
43
44# the main benders loop.
45for i in xsequence(max_iterations):
46
47   print "\nIteration=",i
48
49   # solve the subproblems.
50   solve_all_instances(solver_manager, solver, sub_insts, suffixes=['rc', 'dual'])
51   for instance in sub_insts:
52      print "Profit for scenario="+instance.name+" is "+str(instance.Exp_Stage2_Profit())     
53   print ""
54         
55
56   # if not converged, add store the pricing information from the sub-problem solutions in the master.
57   mstr_inst.CUTS.add(i)
58   for s in xsequence(len(sub_insts)):
59      inst = sub_insts[s-1]
60      for t in mstr_inst.TWOPLUSWEEKS:
61         mstr_inst.time_price[t, s, i] = inst.Time[t].dual
62      for p in mstr_inst.PROD:
63         mstr_inst.bal2_price[p, s, i] = inst.Balance2[p].dual
64      for p in mstr_inst.PROD:
65         for t in mstr_inst.TWOPLUSWEEKS:
66            mstr_inst.sell_lim_price[p, t, s, i] = inst.Sell[p, t].get_suffix_value("Urc")
67
68   # add the master cut.
69   cut = sum([mstr_inst.time_price[t, s, i] * mstr_inst.avail[t] for t in mstr_inst.TWOPLUSWEEKS for s in mstr_inst.SCEN]) + \
70         sum([mstr_inst.bal2_price[p, s, i] * (-mstr_inst.Inv1[p]) for p in mstr_inst.PROD for s in mstr_inst.SCEN]) + \
71         sum([mstr_inst.sell_lim_price[p, t, s, i] * mstr_inst.market[p, t] for p in mstr_inst.PROD for t in mstr_inst.TWOPLUSWEEKS for s in mstr_inst.SCEN]) - \
72         mstr_inst.Min_Stage2_Profit
73   mstr_inst.Cut_Defn.add(i, (0.0, cut, None))
74     
75   # compute expected second-stage profit
76   Exp_Stage2_Profit = sum([sub_insts[s-1].Exp_Stage2_Profit() for s in xsequence(mstr_inst.NUMSCEN())])
77   print "Expected Stage2 Profit=",Exp_Stage2_Profit
78   print ""
79
80   newGAP = mstr_inst.Min_Stage2_Profit.value - Exp_Stage2_Profit
81   print "New gap=",newGAP,"\n"
82
83   if newGAP > 0.00001:
84      GAP = min(GAP, newGAP)
85   else:
86      break
87
88   # re-solve the master and update the subproblem inv1 values.
89   mstr_inst.preprocess() # preprocess to account for the added cuts!
90   solve_all_instances(solver_manager, solver, [mstr_inst])
91
92   print "Master expected profit="+str(mstr_inst.Expected_Profit())
93
94   for instance in sub_insts:
95      for p in mstr_inst.PROD:
96         # the master inventory values might be slightly
97         # less than 0 (within tolerance); threshold here.
98         instance.inv1[p] = max(mstr_inst.Inv1[p](),0.0)
99      instance.preprocess() # re-generate constraints based on new parameter data
100
101print "Benders converged!"
102
103print "\nConverged master solution values:"
104for p in mstr_inst.PROD:
105   print "Make1["+p+"]="+str(mstr_inst.Make1[p]())
106   print "Sell1["+p+"]="+str(mstr_inst.Sell1[p]())   
107   print "Inv1["+p+"]="+str(mstr_inst.Inv1[p]())
Note: See TracBrowser for help on using the repository browser.