benders_2stage.py
Go to the documentation of this file.
18
19from __future__ import print_function
20from gams import *
21import sys
22
24 return '''
25Sets
26 i factories /f1*f3/
27 j distribution centers /d1*d5/
28
29Parameter
30 capacity(i) unit capacity at factories
31 /f1 500, f2 450, f3 650/
32 demand(j) unit demand at distribution centers
33 /d1 160, d2 120, d3 270, d4 325, d5 700 /
34 prodcost unit production cost /14/
35 price sales price /24/
36 wastecost cost of removal of overstocked products /4/
37
38Table transcost(i,j) unit transportation cost
39 d1 d2 d3 d4 d5
40 f1 2.49 5.21 3.76 4.85 2.07
41 f2 1.46 2.54 1.83 1.86 4.76
42 f3 3.26 3.08 2.60 3.76 4.45;
43
44$ifthen not set useBig
45Set
46 s scenarios /lo,mid,hi/
47
48table ScenarioData(s,*) possible outcomes for demand plus probabilities
49 d1 d2 d3 d4 d5 prob
50lo 150 100 250 300 600 0.25
51mid 160 120 270 325 700 0.50
52hi 170 135 300 350 800 0.25;
53$else
54$if not set nrScen $set nrScen 10
55Set s scenarios /s1*s%nrScen%/;
56parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;
57option seed=1234;
58ScenarioData(s,'prob') = 1/card(s);
59ScenarioData(s,j) = demand(j)*uniform(0.6,1.4);
60$endif
61'''
62
64 return '''
65Sets
66 i factories
67 j distribution centers
68
69Parameter
70 capacity(i) unit capacity at factories
71 prodcost unit production cost
72 transcost(i,j) unit transportation cost
73
74$if not set datain $abort 'datain not set'
75$gdxin %datain%
76$load i j capacity prodcost transcost
77
78* Benders master problem
79$if not set maxiter $set maxiter 25
80Set
81 iter max Benders iterations /1*%maxiter%/
82
83Parameter
84 cutconst(iter) constants in optimality cuts
85 cutcoeff(iter,j) coefficients in optimality cuts
86
87Variables
88 ship(i,j) shipments
89 product(i) production
90 received(j) quantity sent to market
91 zmaster objective variable of master problem
92 theta future profit
93Positive Variables ship;
94
95Equations
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));
108product.up(i) = capacity(i);
109
110optcut(iter).. theta =l= cutconst(iter) +
111 sum(j, cutcoeff(iter,j)*received(j));
112
113model masterproblem /all/;
114
115* Initialize cut to be non-binding
116cutconst(iter) = 1e15;
117cutcoeff(iter,j) = eps;
118'''
119
121 return '''
122Sets
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
136* Benders' subproblem
137
138Variables
139 sales(j) sales (actually sold)
140 waste(j) overstocked products
141 zsub objective variable of sub problem
142Positive variables sales, waste
143
144Equations
145 subobj subproblem objective function
146 selling(j) part of received is sold
147 market(j) upperbound on sales
148;
149
150subobj..
151 zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
152
153selling(j).. sales(j) + waste(j) =e= received(j);
154
155market(j).. sales(j) =l= demand(j);
156
157model subproblem /subobj,selling,market/;
158
159* Initialize received
160received(j) = demand(j);
161'''
162
163
164if __name__ == "__main__":
165 if len(sys.argv) > 1:
166 ws = GamsWorkspace(system_directory = sys.argv[1])
167 else:
168 ws = GamsWorkspace()
169
170 data = ws.add_job_from_string(get_data_text())
171
172 opt_data = ws.add_options()
173 opt_data.defines["useBig"] = "1"
174 opt_data.defines["nrScen"] = "100"
175 data.run(opt_data)
176 del opt_data
177
178 scenario_data = data.out_db["ScenarioData"]
179 opt = ws.add_options()
180 opt.defines["datain"] = data.out_db.name
181 maxiter = 40
182 opt.defines["maxiter"] = str(maxiter)
183 opt.all_model_types = "cplexd"
184
185 cp_master = ws.add_checkpoint()
186 cp_sub = ws.add_checkpoint()
187
188 master = ws.add_job_from_string(get_master_text())
189 master.run(opt, cp_master, databases=data.out_db)
190
191 masteri = cp_master.add_modelinstance()
192 cutconst = masteri.sync_db.add_parameter("cutconst", 1, "Benders optimality cut constant")
193 cutcoeff = masteri.sync_db.add_parameter("cutcoeff", 2, "Benders optimality coefficients")
194 theta = masteri.sync_db.add_variable("theta", 0, VarType.Free, "Future profit function variable")
195 theta_fix = masteri.sync_db.add_parameter("thetaFix", 0, "")
196 masteri.instantiate("masterproblem max zmaster using lp", [GamsModifier(cutconst), GamsModifier(cutcoeff), GamsModifier(theta, UpdateAction.Fixed, theta_fix)], opt)
197
198 sub = ws.add_job_from_string(get_sub_text())
199 sub.run(opt, cp_sub, databases=data.out_db)
200 subi = cp_sub.add_modelinstance()
201 received = subi.sync_db.add_parameter("received", 1, "units received from master")
202 demand = subi.sync_db.add_parameter("demand", 1, "stochastic demand")
203 subi.instantiate("subproblem max zsub using lp", [GamsModifier(received), GamsModifier(demand)], opt)
204
205 del opt
206
207 lowerbound = float('-inf')
208 upperbound = float('inf')
209 objmaster = float('inf')
210 iter = 1
211
212 while True:
213 print("Iteration: " + str(iter))
214
215 # Solve master
216 if 1 == iter: # fix theta for first iteration
217 theta_fix.add_record().value = 0
218 else:
219 theta_fix.clear()
220
221 masteri.solve(SymbolUpdateType.BaseCase)
222 print(" Master " + str(masteri.model_status) + " : obj=" + str(masteri.sync_db.get_variable("zmaster").first_record().level))
223 if 1 < iter:
224 upperbound = masteri.sync_db.get_variable("zmaster").first_record().level
225 objmaster = masteri.sync_db.get_variable("zmaster").first_record().level - theta.first_record().level
226
227 # Set received from master
228 received.clear()
229 for r in masteri.sync_db["received"]:
230 received.add_record(r.keys).value = r.level
231 cutcoeff.add_record((str(iter), r.key(0)))
232
233 cutconst.add_record(str(iter))
234 objsub = 0.0
235 for s in data.out_db["s"]:
236 demand.clear()
237 for j in data.out_db["j"]:
238 demand.add_record(j.keys).value = scenario_data.find_record((s.key(0), j.key(0))).value
239
240 subi.solve(SymbolUpdateType.BaseCase)
241 print(" Sub " + str(subi.model_status) + " : obj=" + str(subi.sync_db["zsub"].first_record().level))
242
243 probability = scenario_data.find_record((s.key(0), "prob")).value
244 objsub += probability * subi.sync_db["zsub"].first_record().level
245 for j in data.out_db["j"]:
246 cutconst.find_record(str(iter)).value += probability * subi.sync_db["market"].find_record(j.keys).marginal * demand.find_record(j.keys).value
247 cutcoeff.find_record((str(iter), j.key(0))).value += probability * subi.sync_db["selling"].find_record(j.keys).marginal
248
249 lowerbound = max(lowerbound, objmaster + objsub)
250 iter += 1
251 if iter == maxiter + 1:
252 raise Exception("Benders out of iterations")
253
254 print(" lowerbound: " + str(lowerbound) + " upperbound: " + str(upperbound) + " objmaster: " + str(objmaster))
255 if not ((upperbound - lowerbound) >= 0.001 * (1 + abs(upperbound))):
256 break;
257
258 del masteri
259 del subi