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