28 def scen_solve(i, subi, cutconst, cutcoeff, dem_queue, objsub, queue_lock, io_lock):
34 dem_dict = dem_queue.get()
37 subi[i].sync_db[
"demand"].clear()
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)
44 print(
" Sub " + str(subi[i].model_status) +
" : obj=" + str(subi[i].sync_db[
"zsub"].first_record().level))
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
57 j distribution centers /d1*d5/
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/
68 Table transcost(i,j) unit transportation cost
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;
74 $ifthen not set useBig
76 s scenarios /lo,mid,hi/
78 table ScenarioData(s,*) possible outcomes for demand plus probabilities
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;
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;
88 ScenarioData(s,'prob') = 1/card(s);
89 ScenarioData(s,j) = demand(j)*uniform(0.6,1.4);
97 j distribution centers
100 capacity(i) unit capacity at factories
101 prodcost unit production cost
102 transcost(i,j) unit transportation cost
104 $if not set datain $abort 'datain not set'
106 $load i j capacity prodcost transcost
108 * Benders master problem
109 $if not set maxiter $set maxiter 25
111 iter max Benders iterations /1*%maxiter%/
114 cutconst(iter) constants in optimality cuts
115 cutcoeff(iter,j) coefficients in optimality cuts
119 product(i) production
120 received(j) quantity sent to market
121 zmaster objective variable of master problem
123 Positive Variables ship;
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;
132 zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j))
133 - sum(i,prodcost*product(i));
135 receive(j).. received(j) =e= sum(i, ship(i,j));
137 production(i).. product(i) =e= sum(j, ship(i,j));
138 product.up(i) = capacity(i);
140 optcut(iter).. theta =l= cutconst(iter) +
141 sum(j, cutcoeff(iter,j)*received(j));
143 model masterproblem /all/;
145 * Initialize cut to be non-binding
146 cutconst(iter) = 1e15;
147 cutcoeff(iter,j) = eps;
154 j distribution centers
157 demand(j) unit demand at distribution centers
159 wastecost cost of removal of overstocked products
160 received(j) first stage decision units received
162 $if not set datain $abort 'datain not set'
164 $load i j demand price wastecost
166 * Benders' subproblem
169 sales(j) sales (actually sold)
170 waste(j) overstocked products
171 zsub objective variable of sub problem
172 Positive variables sales, waste
175 subobj subproblem objective function
176 selling(j) part of received is sold
177 market(j) upperbound on sales
181 zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
183 selling(j).. sales(j) + waste(j) =e= received(j);
185 market(j).. sales(j) =l= demand(j);
187 model subproblem /subobj,selling,market/;
189 * Initialize received
190 received(j) = demand(j); '''
193 if __name__ ==
"__main__":
194 if len(sys.argv) > 1:
195 ws = GamsWorkspace(system_directory = sys.argv[1])
201 opt_data = ws.add_options()
202 opt_data.defines[
"useBig"] =
"1"
203 opt_data.defines[
"nrScen"] =
"100"
207 scenario_data = data.out_db.get_parameter(
"ScenarioData")
209 opt = ws.add_options()
210 opt.defines[
"datain"] = data.out_db.name
212 opt.defines[
"maxiter"] = str(maxiter)
213 opt.all_model_types =
"cplex"
215 cp_master = ws.add_checkpoint()
216 cp_sub = ws.add_checkpoint()
219 master.run(opt, cp_master, databases=data.out_db)
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)
228 ws.add_job_from_string(
get_sub_text()).run(opt, cp_sub, databases=data.out_db)
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)
237 for i
in range(1,num_threads):
238 subi.append(subi[0].copy_modelinstance())
241 lowerbound = float(
'-inf')
242 upperbound = float(
'inf')
243 objmaster = float(
'inf')
247 print(
"Iteration: " + str(it))
251 theta_fix.add_record().value = 0
255 masteri.solve(SymbolUpdateType.BaseCase)
256 print(
" Master " + str(masteri.model_status) +
" : obj=" + str(masteri.sync_db[
"zmaster"].first_record().level))
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"]:
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))
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
273 cutconst.add_record(str(it))
276 queue_lock = threading.Lock()
277 io_lock = threading.Lock()
282 for i
in range(num_threads):
286 for j
in data.out_db[
"j"]:
287 coef[i][j.key(0)] = 0.0
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))
294 for i
in range(num_threads):
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)
306 if it == maxiter + 1:
307 raise Exception(
"Benders out of iterations")
309 print(
" lowerbound: " + str(lowerbound) +
" upperbound: " + str(upperbound) +
" objmaster: " + str(objmaster))
311 if not ((upperbound - lowerbound) >= 0.001 * (1 + abs(upperbound))):