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