benders_2stage.py
Go to the documentation of this file.
1
17 
18 from gams import *
19 import sys
20 
21 def get_data_text():
22  return '''
23 Sets
24  i factories /f1*f3/
25  j distribution centers /d1*d5/
26 
27 Parameter
28  capacity(i) unit capacity at factories
29  /f1 500, f2 450, f3 650/
30  demand(j) unit demand at distribution centers
31  /d1 160, d2 120, d3 270, d4 325, d5 700 /
32  prodcost unit production cost /14/
33  price sales price /24/
34  wastecost cost of removal of overstocked products /4/
35 
36 Table transcost(i,j) unit transportation cost
37  d1 d2 d3 d4 d5
38  f1 2.49 5.21 3.76 4.85 2.07
39  f2 1.46 2.54 1.83 1.86 4.76
40  f3 3.26 3.08 2.60 3.76 4.45;
41 
42 $ifthen not set useBig
43 Set
44  s scenarios /lo,mid,hi/
45 
46 table ScenarioData(s,*) possible outcomes for demand plus probabilities
47  d1 d2 d3 d4 d5 prob
48 lo 150 100 250 300 600 0.25
49 mid 160 120 270 325 700 0.50
50 hi 170 135 300 350 800 0.25;
51 $else
52 $if not set nrScen $set nrScen 10
53 Set s scenarios /s1*s%nrScen%/;
54 parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;
55 option seed=1234;
56 ScenarioData(s,'prob') = 1/card(s);
57 ScenarioData(s,j) = demand(j)*uniform(0.6,1.4);
58 $endif
59 '''
60 
61 def get_master_text():
62  return '''
63 Sets
64  i factories
65  j distribution centers
66 
67 Parameter
68  capacity(i) unit capacity at factories
69  prodcost unit production cost
70  transcost(i,j) unit transportation cost
71 
72 $if not set datain $abort 'datain not set'
73 $gdxin %datain%
74 $load i j capacity prodcost transcost
75 
76 * Benders master problem
77 $if not set maxiter $set maxiter 25
78 Set
79  iter max Benders iterations /1*%maxiter%/
80 
81 Parameter
82  cutconst(iter) constants in optimality cuts
83  cutcoeff(iter,j) coefficients in optimality cuts
84 
85 Variables
86  ship(i,j) shipments
87  product(i) production
88  received(j) quantity sent to market
89  zmaster objective variable of master problem
90  theta future profit
91 Positive Variables ship;
92 
93 Equations
94  masterobj master objective function
95  production(i) calculate production in each factory
96  receive(j) calculate quantity to be send to markets
97  optcut(iter) Benders optimality cuts;
98 
99 masterobj..
100  zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j))
101  - sum(i,prodcost*product(i));
102 
103 receive(j).. received(j) =e= sum(i, ship(i,j));
104 
105 production(i).. product(i) =e= sum(j, ship(i,j));
106 product.up(i) = capacity(i);
107 
108 optcut(iter).. theta =l= cutconst(iter) +
109  sum(j, cutcoeff(iter,j)*received(j));
110 
111 model masterproblem /all/;
112 
113 * Initialize cut to be non-binding
114 cutconst(iter) = 1e15;
115 cutcoeff(iter,j) = eps;
116 '''
117 
118 def get_sub_text():
119  return '''
120 Sets
121  i factories
122  j distribution centers
123 
124 Parameter
125  demand(j) unit demand at distribution centers
126  price sales price
127  wastecost cost of removal of overstocked products
128  received(j) first stage decision units received
129 
130 $if not set datain $abort 'datain not set'
131 $gdxin %datain%
132 $load i j demand price wastecost
133 
134 * Benders' subproblem
135 
136 Variables
137  sales(j) sales (actually sold)
138  waste(j) overstocked products
139  zsub objective variable of sub problem
140 Positive variables sales, waste
141 
142 Equations
143  subobj subproblem objective function
144  selling(j) part of received is sold
145  market(j) upperbound on sales
146 ;
147 
148 subobj..
149  zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
150 
151 selling(j).. sales(j) + waste(j) =e= received(j);
152 
153 market(j).. sales(j) =l= demand(j);
154 
155 model subproblem /subobj,selling,market/;
156 
157 * Initialize received
158 received(j) = demand(j);
159 '''
160 
161 
162 if __name__ == "__main__":
163  if len(sys.argv) > 1:
164  ws = GamsWorkspace(system_directory = sys.argv[1])
165  else:
166  ws = GamsWorkspace()
167 
168  data = ws.add_job_from_string(get_data_text())
169 
170  opt_data = ws.add_options()
171  opt_data.defines["useBig"] = "1"
172  opt_data.defines["nrScen"] = "100"
173  data.run(opt_data)
174  del opt_data
175 
176  scenario_data = data.out_db["ScenarioData"]
177  opt = ws.add_options()
178  opt.defines["datain"] = data.out_db.name
179  maxiter = 40
180  opt.defines["maxiter"] = str(maxiter)
181  opt.all_model_types = "cplex"
182 
183  cp_master = ws.add_checkpoint()
184  cp_sub = ws.add_checkpoint()
185 
186  master = ws.add_job_from_string(get_master_text())
187  master.run(opt, cp_master, databases=data.out_db)
188 
189  masteri = cp_master.add_modelinstance()
190  cutconst = masteri.sync_db.add_parameter("cutconst", 1, "Benders optimality cut constant")
191  cutcoeff = masteri.sync_db.add_parameter("cutcoeff", 2, "Benders optimality coefficients")
192  theta = masteri.sync_db.add_variable("theta", 0, VarType.Free, "Future profit function variable")
193  theta_fix = masteri.sync_db.add_parameter("thetaFix", 0, "")
194  masteri.instantiate("masterproblem max zmaster using lp", [GamsModifier(cutconst), GamsModifier(cutcoeff), GamsModifier(theta, UpdateAction.Fixed, theta_fix)], opt)
195 
196  sub = ws.add_job_from_string(get_sub_text())
197  sub.run(opt, cp_sub, databases=data.out_db)
198  subi = cp_sub.add_modelinstance()
199  received = subi.sync_db.add_parameter("received", 1, "units received from master")
200  demand = subi.sync_db.add_parameter("demand", 1, "stochastic demand")
201  subi.instantiate("subproblem max zsub using lp", [GamsModifier(received), GamsModifier(demand)], opt)
202 
203  del opt
204 
205  lowerbound = float('-inf')
206  upperbound = float('inf')
207  objmaster = float('inf')
208  iter = 1
209 
210  while True:
211  print("Iteration: " + str(iter))
212 
213  # Solve master
214  if 1 == iter: # fix theta for first iteration
215  theta_fix.add_record().value = 0
216  else:
217  theta_fix.clear()
218 
219  masteri.solve(SymbolUpdateType.BaseCase)
220  print(" Master " + str(masteri.model_status) + " : obj=" + str(masteri.sync_db.get_variable("zmaster").first_record().level))
221  if 1 < iter:
222  upperbound = masteri.sync_db.get_variable("zmaster").first_record().level
223  objmaster = masteri.sync_db.get_variable("zmaster").first_record().level - theta.first_record().level
224 
225  # Set received from master
226  received.clear()
227  for r in masteri.sync_db["received"]:
228  received.add_record(r.keys).value = r.level
229  cutcoeff.add_record((str(iter), r.key(0)))
230 
231  cutconst.add_record(str(iter))
232  objsub = 0.0
233  for s in data.out_db["s"]:
234  demand.clear()
235  for j in data.out_db["j"]:
236  demand.add_record(j.keys).value = scenario_data.find_record((s.key(0), j.key(0))).value
237 
238  subi.solve(SymbolUpdateType.BaseCase)
239  print(" Sub " + str(subi.model_status) + " : obj=" + str(subi.sync_db["zsub"].first_record().level))
240 
241  probability = scenario_data.find_record((s.key(0), "prob")).value
242  objsub += probability * subi.sync_db["zsub"].first_record().level
243  for j in data.out_db["j"]:
244  cutconst.find_record(str(iter)).value += probability * subi.sync_db["market"].find_record(j.keys).marginal * demand.find_record(j.keys).value
245  cutcoeff.find_record((str(iter), j.key(0))).value += probability * subi.sync_db["selling"].find_record(j.keys).marginal
246 
247  lowerbound = max(lowerbound, objmaster + objsub)
248  iter += 1
249  if iter == maxiter + 1:
250  raise Exception("Benders out of iterations")
251 
252  print(" lowerbound: " + str(lowerbound) + " upperbound: " + str(upperbound) + " objmaster: " + str(objmaster))
253  if not ((upperbound - lowerbound) >= 0.001 * (1 + abs(upperbound))):
254  break;
255 
256  del masteri
257  del subi
benders_2stage.get_master_text
def get_master_text()
Definition: benders_2stage.py:61
benders_2stage.get_data_text
def get_data_text()
Definition: benders_2stage.py:21
benders_2stage.get_sub_text
def get_sub_text()
Definition: benders_2stage.py:118