20from queue
import Queue
22from threading
import Lock, Thread
23from gams
import GamsModifier, GamsWorkspace, SymbolUpdateType, UpdateAction, VarType
27 i 'factories' / f1*f3 /
28 j
'distribution centers' / d1*d5 /;
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 /;
39Table transcost(i,j)
'unit transportation cost'
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;
46 Set s
'scenarios' / lo, mid, hi /;
48 Table ScenarioData(s,*)
'possible outcomes for demand plus probabilities'
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;
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';
58 ScenarioData(s,
'prob') = 1/card(s);
59 ScenarioData(s,j) =
demand(j)*uniform(0.6,1.4);
63GAMS_MASTER_MODEL = """
66 j
'distribution centers';
69 capacity(i)
'unit capacity at factories'
70 prodcost
'unit production cost'
71 transcost(i,j)
'unit transportation cost';
73$
if not set datain $abort
'datain not set'
75$load i j capacity prodcost transcost
78* Benders master problem
79$
if not set max_iter $set max_iter 25
80Set iter
'max Benders iterations' / 1*%max_iter% /;
83 cutconst(iter)
'constants in optimality cuts'
84 cutcoeff(iter,j)
'coefficients in optimality cuts';
88 product(i)
'production'
89 received(j)
'quantity sent to market'
90 zmaster
'objective variable of master problem'
91 theta
'future profit';
93Positive Variable ship;
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';
102 zmaster =e= theta - sum((i,j), transcost(i,j)*ship(i,j))
103 - sum(i, prodcost*product(i));
105receive(j)..
received(j) =e= sum(i, ship(i,j));
107production(i).. product(i) =e= sum(j, ship(i,j));
112product.up(i) = capacity(i);
114Model masterproblem / all /;
116* Initialize cut to be non-binding
124 j
'distribution centers';
127 demand(j)
'unit demand at distribution centers'
129 wastecost
'cost of removal of overstocked products'
130 received(j)
'first stage decision units received';
132$
if not set datain $abort
'datain not set'
134$load i j demand price wastecost
139 sales(j) 'sales (actually sold)'
140 waste(j)
'overstocked products'
141 zsub
'objective variable of sub problem';
143Positive Variable sales, waste;
146 subobj
'subproblem objective function'
147 selling(j)
'part of received is sold'
148 market(j)
'upper_bound on sales';
150subobj.. zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
152selling(j).. sales(j) + waste(j) =e=
received(j);
154market(j).. sales(j) =l=
demand(j);
156Model subproblem / subobj, selling, market /;
163def scen_solve(i, mi_sub, cutconst, cutcoeff, dem_queue, obj_sub, queue_lock, io_lock):
166 if dem_queue.empty():
169 dem_dict = dem_queue.get()
172 mi_sub[i].sync_db[
"demand"].clear()
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)
180 f
" Sub {mi_sub[i].model_status} : obj={mi_sub[i].sync_db['zsub'].first_record().level}"
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()):
188 probability * mi_sub[i].sync_db[
"market"].find_record(k).marginal * v
191 probability * mi_sub[i].sync_db[
"selling"].find_record(k).marginal
195if __name__ ==
"__main__":
196 sys_dir = sys.argv[1]
if len(sys.argv) > 1
else None
197 ws = GamsWorkspace(system_directory=sys_dir)
199 data = ws.add_job_from_string(GAMS_DATA)
201 opt_data = ws.add_options()
202 opt_data.defines[
"useBig"] =
"1"
203 opt_data.defines[
"nrScen"] =
"100"
206 scenario_data = data.out_db[
"ScenarioData"]
207 opt = ws.add_options()
208 opt.defines[
"datain"] = data.out_db.name
210 opt.defines[
"max_iter"] = str(max_iter)
211 opt.all_model_types =
"cplex"
213 cp_master = ws.add_checkpoint()
214 cp_sub = ws.add_checkpoint()
216 master = ws.add_job_from_string(GAMS_MASTER_MODEL)
217 master.run(opt, cp_master, databases=data.out_db)
219 mi_master = cp_master.add_modelinstance()
220 cutconst = mi_master.sync_db.add_parameter(
221 "cutconst", 1,
"Benders optimality cut constant"
223 cutcoeff = mi_master.sync_db.add_parameter(
224 "cutcoeff", 2,
"Benders optimality coefficients"
226 theta = mi_master.sync_db.add_variable(
227 "theta", 0, VarType.Free,
"Future profit function variable"
229 theta_fix = mi_master.sync_db.add_parameter(
"thetaFix", 0)
230 mi_master.instantiate(
231 "masterproblem max zmaster using lp",
233 GamsModifier(cutconst),
234 GamsModifier(cutcoeff),
235 GamsModifier(theta, UpdateAction.Fixed, theta_fix),
240 sub = ws.add_job_from_string(GAMS_SUB_MODEL)
241 sub.run(opt, cp_sub, databases=data.out_db)
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"
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)],
256 mi_sub += [mi_sub[0].copy_modelinstance()
for i
in range(num_threads - 1)]
258 lower_bound = float(
"-inf")
259 upper_bound = float(
"inf")
260 obj_master = float(
"inf")
262 theta_fix.add_record().value = 0
265 print(f
"Iteration: {it}")
271 mi_master.solve(SymbolUpdateType.BaseCase)
273 f
" Master {mi_master.model_status} : obj={mi_master.sync_db.get_variable('zmaster').first_record().level}"
276 upper_bound = mi_master.sync_db[
"zmaster"].first_record().level
278 mi_master.sync_db[
"zmaster"].first_record().level
279 - theta.first_record().level
281 for s
in data.out_db[
"s"]:
283 j.key(0): scenario_data.find_record([s.key(0), j.key(0)]).value
284 for j
in data.out_db[
"j"]
289 scenario_data.find_record([s.key(0),
"prob"]).value,
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
301 cutconst.add_record(str(it))
305 obj_sub = [0.0] * num_threads
306 cons = [0.0] * num_threads
308 i: {j.key(0): 0.0
for j
in data.out_db[
"j"]}
for i
in range(num_threads)
313 for i
in range(num_threads):
316 args=(i, mi_sub, cons, coef, dem_queue, obj_sub, queue_lock, io_lock),
319 for i
in range(num_threads):
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)]
329 lower_bound = max(lower_bound, obj_master + obj_sub_sum)
332 raise Exception(
"Benders out of iterations")
335 f
" lower bound: {lower_bound} upper bound: {upper_bound} obj master: {obj_master}"
337 if upper_bound - lower_bound < 0.001 * (1 + abs(upper_bound)):