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