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 ws = GamsWorkspace(system_directory=sys_dir)
198
199 data = ws.add_job_from_string(GAMS_DATA)
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
206 scenario_data = data.out_db["ScenarioData"]
207 opt = ws.add_options()
208 opt.defines["datain"] = data.out_db.name
209 max_iter = 40
210 opt.defines["max_iter"] = str(max_iter)
211 opt.all_model_types = "cplex"
212
213 cp_master = ws.add_checkpoint()
214 cp_sub = ws.add_checkpoint()
215
216 master = ws.add_job_from_string(GAMS_MASTER_MODEL)
217 master.run(opt, cp_master, databases=data.out_db)
218
219 mi_master = cp_master.add_modelinstance()
220 cutconst = mi_master.sync_db.add_parameter(
221 "cutconst", 1, "Benders optimality cut constant"
222 )
223 cutcoeff = mi_master.sync_db.add_parameter(
224 "cutcoeff", 2, "Benders optimality coefficients"
225 )
226 theta = mi_master.sync_db.add_variable(
227 "theta", 0, VarType.Free, "Future profit function variable"
228 )
229 theta_fix = mi_master.sync_db.add_parameter("thetaFix", 0)
230 mi_master.instantiate(
231 "masterproblem max zmaster using lp",
232 [
233 GamsModifier(cutconst),
234 GamsModifier(cutcoeff),
235 GamsModifier(theta, UpdateAction.Fixed, theta_fix),
236 ],
237 opt,
238 )
239
240 sub = ws.add_job_from_string(GAMS_SUB_MODEL)
241 sub.run(opt, cp_sub, databases=data.out_db)
242 num_threads = 2
243 mi_sub = []
244 dem_queue = Queue()
245 mi_sub.append(cp_sub.add_modelinstance())
246 received = mi_sub[0].sync_db.add_parameter(
247 "received", 1, "units received from first stage solution"
248 )
249 demand = mi_sub[0].sync_db.add_parameter("demand", 1, "stochastic demand")
250 mi_sub[0].instantiate(
251 "subproblem max zsub using lp",
252 [GamsModifier(received), GamsModifier(demand)],
253 opt,
254 )
255
256 mi_sub += [mi_sub[0].copy_modelinstance() for i in range(num_threads - 1)]
257
258 lower_bound = float("-inf")
259 upper_bound = float("inf")
260 obj_master = float("inf")
261 it = 1
262 theta_fix.add_record().value = 0 # fix theta for first iteration
263
264 while True:
265 print(f"Iteration: {it}")
266
267 # solve master
268 if it > 1:
269 theta_fix.clear()
270
271 mi_master.solve(SymbolUpdateType.BaseCase)
272 print(
273 f" Master {mi_master.model_status} : obj={mi_master.sync_db.get_variable('zmaster').first_record().level}"
274 )
275 if it > 1:
276 upper_bound = mi_master.sync_db["zmaster"].first_record().level
277 obj_master = (
278 mi_master.sync_db["zmaster"].first_record().level
279 - theta.first_record().level
280 )
281 for s in data.out_db["s"]:
282 dem_dict = {
283 j.key(0): scenario_data.find_record([s.key(0), j.key(0)]).value
284 for j in data.out_db["j"]
285 }
286 dem_queue.put(
287 (
288 s.key(0),
289 scenario_data.find_record([s.key(0), "prob"]).value,
290 dem_dict,
291 )
292 )
293
294 for i in range(num_threads):
295 mi_sub[i].sync_db["received"].clear()
296 for r in mi_master.sync_db["received"]:
297 cutcoeff.add_record([str(it), r.key(0)])
298 for i in range(num_threads):
299 mi_sub[i].sync_db["received"].add_record(r.keys).value = r.level
300
301 cutconst.add_record(str(it))
302
303 queue_lock = Lock()
304 io_lock = Lock()
305 obj_sub = [0.0] * num_threads
306 cons = [0.0] * num_threads
307 coef = {
308 i: {j.key(0): 0.0 for j in data.out_db["j"]} for i in range(num_threads)
309 }
310
311 # solve multiple model instances in parallel
312 threads = {}
313 for i in range(num_threads):
314 threads[i] = Thread(
315 target=scen_solve,
316 args=(i, mi_sub, cons, coef, dem_queue, obj_sub, queue_lock, io_lock),
317 )
318 threads[i].start()
319 for i in range(num_threads):
320 threads[i].join()
321
322 obj_sub_sum = 0.0
323 for i in range(num_threads):
324 obj_sub_sum += obj_sub[i]
325 cutconst.find_record(str(it)).value += cons[i]
326 for j in data.out_db["j"]:
327 cutcoeff.find_record([str(it), j.key(0)]).value += coef[i][j.key(0)]
328
329 lower_bound = max(lower_bound, obj_master + obj_sub_sum)
330 it += 1
331 if it > max_iter:
332 raise Exception("Benders out of iterations")
333
334 print(
335 f" lower bound: {lower_bound} upper bound: {upper_bound} obj master: {obj_master}"
336 )
337 if upper_bound - lower_bound < 0.001 * (1 + abs(upper_bound)):
338 break
benders_2stage_mt.received
list received
Definition: benders_2stage_mt.py:246
benders_2stage_mt.demand
list demand
Definition: benders_2stage_mt.py:249
benders_2stage_mt.cutcoeff
GamsWorkspace cutcoeff
Definition: benders_2stage_mt.py:223
benders_2stage_mt.cutconst
GamsWorkspace cutconst
Definition: benders_2stage_mt.py:220