18from gams
import GamsModifier, GamsWorkspace, SymbolUpdateType, UpdateAction, VarType
22 i 'factories' / f1*f3 /
23 j
'distribution centers' / d1*d5 /;
26 capacity(i)
'unit capacity at factories'
27 / f1 500, f2 450, f3 650 /
28 demand(j)
'unit demand at distribution centers'
29 / d1 160, d2 120, d3 270, d4 325, d5 700 /
30 prodcost
'unit production cost' / 14 /
31 price
'sales price' / 24 /
32 wastecost
'cost of removal of overstocked products' / 4 /;
34Table transcost(i,j)
'unit transportation cost'
36 f1 2.49 5.21 3.76 4.85 2.07
37 f2 1.46 2.54 1.83 1.86 4.76
38 f3 3.26 3.08 2.60 3.76 4.45;
41 Set s
'scenarios' / lo, mid, hi /;
43 Table ScenarioData(s,*)
'possible outcomes for demand plus probabilities'
45 lo 150 100 250 300 600 0.25
46 mid 160 120 270 325 700 0.50
47 hi 170 135 300 350 800 0.25;
49$
if not set nrScen $set nrScen 10
50 Set s
'scenarios' / s1*s%nrScen% /;
51 Parameter ScenarioData(s,*)
'possible outcomes for demand plus probabilities';
53 ScenarioData(s,
'prob') = 1/card(s);
54 ScenarioData(s,j) =
demand(j)*uniform(0.6,1.4);
58GAMS_MASTER_MODEL = """
61 j
'distribution centers';
64 capacity(i)
'unit capacity at factories'
65 prodcost
'unit production cost'
66 transcost(i,j)
'unit transportation cost';
68$
if not set datain $abort
'datain not set'
70$load i j capacity prodcost transcost
73* Benders master problem
74$
if not set max_iter $set max_iter 25
75Set iter
'max Benders iterations' / 1*%max_iter% /;
78 cutconst(iter)
'constants in optimality cuts'
79 cutcoeff(iter,j)
'coefficients in optimality cuts';
83 product(i)
'production'
84 received(j)
'quantity sent to market'
85 zmaster
'objective variable of master problem'
86 theta
'future profit';
88Positive Variable ship;
91 masterobj
'master objective function'
92 production(i)
'calculate production in each factory'
93 receive(j)
'calculate quantity to be send to markets'
94 optcut(iter)
'Benders optimality cuts';
97 zmaster =e= theta - sum((i,j), transcost(i,j)*ship(i,j))
98 - sum(i, prodcost*product(i));
100receive(j)..
received(j) =e= sum(i, ship(i,j));
102production(i).. product(i) =e= sum(j, ship(i,j));
107product.up(i) = capacity(i);
109Model masterproblem / all /;
111* Initialize cut to be non-binding
119 j
'distribution centers';
122 demand(j)
'unit demand at distribution centers'
124 wastecost
'cost of removal of overstocked products'
125 received(j)
'first stage decision units received';
127$
if not set datain $abort
'datain not set'
129$load i j demand price wastecost
134 sales(j) 'sales (actually sold)'
135 waste(j)
'overstocked products'
136 zsub
'objective variable of sub problem';
138Positive Variable sales, waste;
141 subobj
'subproblem objective function'
142 selling(j)
'part of received is sold'
143 market(j)
'upper_bound on sales';
145subobj.. zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
147selling(j).. sales(j) + waste(j) =e=
received(j);
149market(j).. sales(j) =l=
demand(j);
151Model subproblem / subobj, selling, market /;
157if __name__ == "__main__":
158 sys_dir = sys.argv[1] if len(sys.argv) > 1
else None
159 ws = GamsWorkspace(system_directory=sys_dir)
161 data = ws.add_job_from_string(GAMS_DATA)
163 opt_data = ws.add_options()
164 opt_data.defines[
"useBig"] =
"1"
165 opt_data.defines[
"nrScen"] =
"100"
168 scenario_data = data.out_db[
"ScenarioData"]
169 opt = ws.add_options()
170 opt.defines[
"datain"] = data.out_db.name
172 opt.defines[
"max_iter"] = str(max_iter)
173 opt.all_model_types =
"cplex"
175 cp_master = ws.add_checkpoint()
176 cp_sub = ws.add_checkpoint()
178 master = ws.add_job_from_string(GAMS_MASTER_MODEL)
179 master.run(opt, cp_master, databases=data.out_db)
181 mi_master = cp_master.add_modelinstance()
182 cutconst = mi_master.sync_db.add_parameter(
183 "cutconst", 1,
"Benders optimality cut constant"
185 cutcoeff = mi_master.sync_db.add_parameter(
186 "cutcoeff", 2,
"Benders optimality coefficients"
188 theta = mi_master.sync_db.add_variable(
189 "theta", 0, VarType.Free,
"Future profit function variable"
191 theta_fix = mi_master.sync_db.add_parameter(
"thetaFix", 0)
192 mi_master.instantiate(
193 "masterproblem max zmaster using lp",
195 GamsModifier(cutconst),
196 GamsModifier(cutcoeff),
197 GamsModifier(theta, UpdateAction.Fixed, theta_fix),
202 sub = ws.add_job_from_string(GAMS_SUB_MODEL)
203 sub.run(opt, cp_sub, databases=data.out_db)
204 mi_sub = cp_sub.add_modelinstance()
205 received = mi_sub.sync_db.add_parameter(
"received", 1,
"units received from master")
206 demand = mi_sub.sync_db.add_parameter(
"demand", 1,
"stochastic demand")
208 "subproblem max zsub using lp",
209 [GamsModifier(received), GamsModifier(demand)],
213 lower_bound = float(
"-inf")
214 upper_bound = float(
"inf")
215 obj_master = float(
"inf")
217 theta_fix.add_record().value = 0
220 print(f
"Iteration: {it}")
226 mi_master.solve(SymbolUpdateType.BaseCase)
228 f
" Master {mi_master.model_status} : obj={mi_master.sync_db.get_variable('zmaster').first_record().level}"
231 upper_bound = mi_master.sync_db.get_variable(
"zmaster").first_record().level
233 mi_master.sync_db[
"zmaster"].first_record().level
234 - theta.first_record().level
239 for r
in mi_master.sync_db[
"received"]:
240 received.add_record(r.keys).value = r.level
241 cutcoeff.add_record((str(it), r.key(0)))
243 cutconst.add_record(str(it))
245 for s
in data.out_db[
"s"]:
247 for j
in data.out_db[
"j"]:
248 demand.add_record(j.keys).value = scenario_data.find_record(
251 mi_sub.solve(SymbolUpdateType.BaseCase)
253 f
" Sub {mi_sub.model_status} : obj={mi_sub.sync_db['zsub'].first_record().level}"
255 probability = scenario_data.find_record((s.key(0),
"prob")).value
256 obj_sub += probability * mi_sub.sync_db[
"zsub"].first_record().level
257 for j
in data.out_db[
"j"]:
258 cutconst.find_record(str(it)).value += (
260 * mi_sub.sync_db[
"market"].find_record(j.keys).marginal
261 * demand.find_record(j.keys).value
263 cutcoeff.find_record((str(it), j.key(0))).value += (
264 probability * mi_sub.sync_db[
"selling"].find_record(j.keys).marginal
267 lower_bound = max(lower_bound, obj_master + obj_sub)
270 raise Exception(
"Benders out of iterations")
273 f
" lower bound: {lower_bound} upper bound: {upper_bound} obj master: {obj_master}"
275 if upper_bound - lower_bound < 0.001 * (1 + abs(upper_bound)):