25 j distribution centers /d1*d5/
28 capacity(i) unit capacity at factories
29 /f1 500, f2 450, f3 650/
30 demand(j) unit demand at distribution centers
31 /d1 160, d2 120, d3 270, d4 325, d5 700 /
32 prodcost unit production cost /14/
33 price sales price /24/
34 wastecost cost of removal of overstocked products /4/
36 Table transcost(i,j) unit transportation cost
38 f1 2.49 5.21 3.76 4.85 2.07
39 f2 1.46 2.54 1.83 1.86 4.76
40 f3 3.26 3.08 2.60 3.76 4.45;
42 $ifthen not set useBig
44 s scenarios /lo,mid,hi/
46 table ScenarioData(s,*) possible outcomes for demand plus probabilities
48 lo 150 100 250 300 600 0.25
49 mid 160 120 270 325 700 0.50
50 hi 170 135 300 350 800 0.25;
52 $if not set nrScen $set nrScen 10
53 Set s scenarios /s1*s%nrScen%/;
54 parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;
56 ScenarioData(s,'prob') = 1/card(s);
57 ScenarioData(s,j) = demand(j)*uniform(0.6,1.4);
65 j distribution centers
68 capacity(i) unit capacity at factories
69 prodcost unit production cost
70 transcost(i,j) unit transportation cost
72 $if not set datain $abort 'datain not set'
74 $load i j capacity prodcost transcost
76 * Benders master problem
77 $if not set maxiter $set maxiter 25
79 iter max Benders iterations /1*%maxiter%/
82 cutconst(iter) constants in optimality cuts
83 cutcoeff(iter,j) coefficients in optimality cuts
88 received(j) quantity sent to market
89 zmaster objective variable of master problem
91 Positive Variables ship;
94 masterobj master objective function
95 production(i) calculate production in each factory
96 receive(j) calculate quantity to be send to markets
97 optcut(iter) Benders optimality cuts;
100 zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j))
101 - sum(i,prodcost*product(i));
103 receive(j).. received(j) =e= sum(i, ship(i,j));
105 production(i).. product(i) =e= sum(j, ship(i,j));
106 product.up(i) = capacity(i);
108 optcut(iter).. theta =l= cutconst(iter) +
109 sum(j, cutcoeff(iter,j)*received(j));
111 model masterproblem /all/;
113 * Initialize cut to be non-binding
114 cutconst(iter) = 1e15;
115 cutcoeff(iter,j) = eps;
122 j distribution centers
125 demand(j) unit demand at distribution centers
127 wastecost cost of removal of overstocked products
128 received(j) first stage decision units received
130 $if not set datain $abort 'datain not set'
132 $load i j demand price wastecost
134 * Benders' subproblem
137 sales(j) sales (actually sold)
138 waste(j) overstocked products
139 zsub objective variable of sub problem
140 Positive variables sales, waste
143 subobj subproblem objective function
144 selling(j) part of received is sold
145 market(j) upperbound on sales
149 zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
151 selling(j).. sales(j) + waste(j) =e= received(j);
153 market(j).. sales(j) =l= demand(j);
155 model subproblem /subobj,selling,market/;
157 * Initialize received
158 received(j) = demand(j);
162 if __name__ ==
"__main__":
163 if len(sys.argv) > 1:
164 ws = GamsWorkspace(system_directory = sys.argv[1])
170 opt_data = ws.add_options()
171 opt_data.defines[
"useBig"] =
"1"
172 opt_data.defines[
"nrScen"] =
"100"
176 scenario_data = data.out_db[
"ScenarioData"]
177 opt = ws.add_options()
178 opt.defines[
"datain"] = data.out_db.name
180 opt.defines[
"maxiter"] = str(maxiter)
181 opt.all_model_types =
"cplex"
183 cp_master = ws.add_checkpoint()
184 cp_sub = ws.add_checkpoint()
187 master.run(opt, cp_master, databases=data.out_db)
189 masteri = cp_master.add_modelinstance()
190 cutconst = masteri.sync_db.add_parameter(
"cutconst", 1,
"Benders optimality cut constant")
191 cutcoeff = masteri.sync_db.add_parameter(
"cutcoeff", 2,
"Benders optimality coefficients")
192 theta = masteri.sync_db.add_variable(
"theta", 0, VarType.Free,
"Future profit function variable")
193 theta_fix = masteri.sync_db.add_parameter(
"thetaFix", 0,
"")
194 masteri.instantiate(
"masterproblem max zmaster using lp", [GamsModifier(cutconst), GamsModifier(cutcoeff), GamsModifier(theta, UpdateAction.Fixed, theta_fix)], opt)
197 sub.run(opt, cp_sub, databases=data.out_db)
198 subi = cp_sub.add_modelinstance()
199 received = subi.sync_db.add_parameter(
"received", 1,
"units received from master")
200 demand = subi.sync_db.add_parameter(
"demand", 1,
"stochastic demand")
201 subi.instantiate(
"subproblem max zsub using lp", [GamsModifier(received), GamsModifier(demand)], opt)
205 lowerbound = float(
'-inf')
206 upperbound = float(
'inf')
207 objmaster = float(
'inf')
211 print(
"Iteration: " + str(iter))
215 theta_fix.add_record().value = 0
219 masteri.solve(SymbolUpdateType.BaseCase)
220 print(
" Master " + str(masteri.model_status) +
" : obj=" + str(masteri.sync_db.get_variable(
"zmaster").first_record().level))
222 upperbound = masteri.sync_db.get_variable(
"zmaster").first_record().level
223 objmaster = masteri.sync_db.get_variable(
"zmaster").first_record().level - theta.first_record().level
227 for r
in masteri.sync_db[
"received"]:
228 received.add_record(r.keys).value = r.level
229 cutcoeff.add_record((str(iter), r.key(0)))
231 cutconst.add_record(str(iter))
233 for s
in data.out_db[
"s"]:
235 for j
in data.out_db[
"j"]:
236 demand.add_record(j.keys).value = scenario_data.find_record((s.key(0), j.key(0))).value
238 subi.solve(SymbolUpdateType.BaseCase)
239 print(
" Sub " + str(subi.model_status) +
" : obj=" + str(subi.sync_db[
"zsub"].first_record().level))
241 probability = scenario_data.find_record((s.key(0),
"prob")).value
242 objsub += probability * subi.sync_db[
"zsub"].first_record().level
243 for j
in data.out_db[
"j"]:
244 cutconst.find_record(str(iter)).value += probability * subi.sync_db[
"market"].find_record(j.keys).marginal * demand.find_record(j.keys).value
245 cutcoeff.find_record((str(iter), j.key(0))).value += probability * subi.sync_db[
"selling"].find_record(j.keys).marginal
247 lowerbound = max(lowerbound, objmaster + objsub)
249 if iter == maxiter + 1:
250 raise Exception(
"Benders out of iterations")
252 print(
" lowerbound: " + str(lowerbound) +
" upperbound: " + str(upperbound) +
" objmaster: " + str(objmaster))
253 if not ((upperbound - lowerbound) >= 0.001 * (1 + abs(upperbound))):