Loading...
Searching...
No Matches
benders_2stage.py
Go to the documentation of this file.
1
16
17import sys
18from gams import GamsModifier, GamsWorkspace, SymbolUpdateType, UpdateAction, VarType
19
20GAMS_DATA = """
21Set
22 i 'factories' / f1*f3 /
23 j 'distribution centers' / d1*d5 /;
24
25Parameter
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 /;
33
34Table transcost(i,j) 'unit transportation cost'
35 d1 d2 d3 d4 d5
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;
39
40$ifThen not set useBig
41 Set s 'scenarios' / lo, mid, hi /;
42
43 Table ScenarioData(s,*) 'possible outcomes for demand plus probabilities'
44 d1 d2 d3 d4 d5 prob
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;
48$else
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';
52 option seed = 1234;
53 ScenarioData(s,'prob') = 1/card(s);
54 ScenarioData(s,j) = demand(j)*uniform(0.6,1.4);
55$endIf
56"""
57
58GAMS_MASTER_MODEL = """
59Set
60 i 'factories'
61 j 'distribution centers';
62
63Parameter
64 capacity(i) 'unit capacity at factories'
65 prodcost 'unit production cost'
66 transcost(i,j) 'unit transportation cost';
67
68$if not set datain $abort 'datain not set'
69$gdxIn %datain%
70$load i j capacity prodcost transcost
71$gdxIn
72
73* Benders master problem
74$if not set max_iter $set max_iter 25
75Set iter 'max Benders iterations' / 1*%max_iter% /;
76
77Parameter
78 cutconst(iter) 'constants in optimality cuts'
79 cutcoeff(iter,j) 'coefficients in optimality cuts';
80
81Variable
82 ship(i,j) 'shipments'
83 product(i) 'production'
84 received(j) 'quantity sent to market'
85 zmaster 'objective variable of master problem'
86 theta 'future profit';
87
88Positive Variable ship;
89
90Equation
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';
95
96masterobj..
97 zmaster =e= theta - sum((i,j), transcost(i,j)*ship(i,j))
98 - sum(i, prodcost*product(i));
99
100receive(j).. received(j) =e= sum(i, ship(i,j));
101
102production(i).. product(i) =e= sum(j, ship(i,j));
103
104optcut(iter)..
105 theta =l= cutconst(iter) + sum(j, cutcoeff(iter,j)*received(j));
106
107product.up(i) = capacity(i);
108
109Model masterproblem / all /;
110
111* Initialize cut to be non-binding
112cutconst(iter) = 1e15;
113cutcoeff(iter,j) = eps;
114"""
115
116GAMS_SUB_MODEL = """
117Set
118 i 'factories'
119 j 'distribution centers';
120
121Parameter
122 demand(j) 'unit demand at distribution centers'
123 price 'sales price'
124 wastecost 'cost of removal of overstocked products'
125 received(j) 'first stage decision units received';
126
127$if not set datain $abort 'datain not set'
128$gdxIn %datain%
129$load i j demand price wastecost
130$gdxIn
131
132* Benders' subproblem
133Variable
134 sales(j) 'sales (actually sold)'
135 waste(j) 'overstocked products'
136 zsub 'objective variable of sub problem';
137
138Positive Variable sales, waste;
139
140Equation
141 subobj 'subproblem objective function'
142 selling(j) 'part of received is sold'
143 market(j) 'upper_bound on sales';
144
145subobj.. zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
146
147selling(j).. sales(j) + waste(j) =e= received(j);
148
149market(j).. sales(j) =l= demand(j);
150
151Model subproblem / subobj, selling, market /;
152
153* Initialize received
154received(j) = demand(j);
155"""
156
157if __name__ == "__main__":
158 sys_dir = sys.argv[1] if len(sys.argv) > 1 else None
159 work_dir = sys.argv[2] if len(sys.argv) > 2 else None
160 ws = GamsWorkspace(system_directory=sys_dir, working_directory=work_dir)
161
162 data = ws.add_job_from_string(GAMS_DATA)
163
164 opt_data = ws.add_options()
165 opt_data.defines["useBig"] = "1"
166 opt_data.defines["nrScen"] = "100"
167 data.run(opt_data)
168
169 scenario_data = data.out_db["ScenarioData"]
170 opt = ws.add_options()
171 opt.defines["datain"] = data.out_db.name
172 max_iter = 40
173 opt.defines["max_iter"] = str(max_iter)
174 opt.all_model_types = "cplex"
175
176 cp_master = ws.add_checkpoint()
177 cp_sub = ws.add_checkpoint()
178
179 master = ws.add_job_from_string(GAMS_MASTER_MODEL)
180 master.run(opt, cp_master, databases=data.out_db)
181
182 mi_master = cp_master.add_modelinstance()
183 cutconst = mi_master.sync_db.add_parameter(
184 "cutconst", 1, "Benders optimality cut constant"
185 )
186 cutcoeff = mi_master.sync_db.add_parameter(
187 "cutcoeff", 2, "Benders optimality coefficients"
188 )
189 theta = mi_master.sync_db.add_variable(
190 "theta", 0, VarType.Free, "Future profit function variable"
191 )
192 theta_fix = mi_master.sync_db.add_parameter("thetaFix", 0)
193 mi_master.instantiate(
194 "masterproblem max zmaster using lp",
195 [
196 GamsModifier(cutconst),
197 GamsModifier(cutcoeff),
198 GamsModifier(theta, UpdateAction.Fixed, theta_fix),
199 ],
200 opt,
201 )
202
203 sub = ws.add_job_from_string(GAMS_SUB_MODEL)
204 sub.run(opt, cp_sub, databases=data.out_db)
205 mi_sub = cp_sub.add_modelinstance()
206 received = mi_sub.sync_db.add_parameter("received", 1, "units received from master")
207 demand = mi_sub.sync_db.add_parameter("demand", 1, "stochastic demand")
208 mi_sub.instantiate(
209 "subproblem max zsub using lp",
210 [GamsModifier(received), GamsModifier(demand)],
211 opt,
212 )
213
214 lower_bound = float("-inf")
215 upper_bound = float("inf")
216 obj_master = float("inf")
217 it = 1
218 theta_fix.add_record().value = 0 # fix theta for first iteration
219
220 while True:
221 print(f"Iteration: {it}")
222
223 # solve master
224 if it > 1:
225 theta_fix.clear()
226
227 mi_master.solve(SymbolUpdateType.BaseCase)
228 print(
229 f" Master {mi_master.model_status} : obj={mi_master.sync_db.get_variable('zmaster').first_record().level}"
230 )
231 if it > 1:
232 upper_bound = mi_master.sync_db.get_variable("zmaster").first_record().level
233 obj_master = (
234 mi_master.sync_db["zmaster"].first_record().level
235 - theta.first_record().level
236 )
237
238 # set received from master
239 received.clear()
240 for r in mi_master.sync_db["received"]:
241 received.add_record(r.keys).value = r.level
242 cutcoeff.add_record((str(it), r.key(0)))
243
244 cutconst.add_record(str(it))
245 obj_sub = 0.0
246 for s in data.out_db["s"]:
247 demand.clear()
248 for j in data.out_db["j"]:
249 demand.add_record(j.keys).value = scenario_data.find_record(
250 (s.key(0), j.key(0))
251 ).value
252 mi_sub.solve(SymbolUpdateType.BaseCase)
253 print(
254 f" Sub {mi_sub.model_status} : obj={mi_sub.sync_db['zsub'].first_record().level}"
255 )
256 probability = scenario_data.find_record((s.key(0), "prob")).value
257 obj_sub += probability * mi_sub.sync_db["zsub"].first_record().level
258 for j in data.out_db["j"]:
259 cutconst.find_record(str(it)).value += (
260 probability
261 * mi_sub.sync_db["market"].find_record(j.keys).marginal
262 * demand.find_record(j.keys).value
263 )
264 cutcoeff.find_record((str(it), j.key(0))).value += (
265 probability * mi_sub.sync_db["selling"].find_record(j.keys).marginal
266 )
267
268 lower_bound = max(lower_bound, obj_master + obj_sub)
269 it += 1
270 if it > max_iter:
271 raise Exception("Benders out of iterations")
272
273 print(
274 f" lower bound: {lower_bound} upper bound: {upper_bound} obj master: {obj_master}"
275 )
276 if upper_bound - lower_bound < 0.001 * (1 + abs(upper_bound)):
277 break