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 ws = GamsWorkspace(system_directory=sys_dir)
160
161 data = ws.add_job_from_string(GAMS_DATA)
162
163 opt_data = ws.add_options()
164 opt_data.defines["useBig"] = "1"
165 opt_data.defines["nrScen"] = "100"
166 data.run(opt_data)
167
168 scenario_data = data.out_db["ScenarioData"]
169 opt = ws.add_options()
170 opt.defines["datain"] = data.out_db.name
171 max_iter = 40
172 opt.defines["max_iter"] = str(max_iter)
173 opt.all_model_types = "cplex"
174
175 cp_master = ws.add_checkpoint()
176 cp_sub = ws.add_checkpoint()
177
178 master = ws.add_job_from_string(GAMS_MASTER_MODEL)
179 master.run(opt, cp_master, databases=data.out_db)
180
181 mi_master = cp_master.add_modelinstance()
182 cutconst = mi_master.sync_db.add_parameter(
183 "cutconst", 1, "Benders optimality cut constant"
184 )
185 cutcoeff = mi_master.sync_db.add_parameter(
186 "cutcoeff", 2, "Benders optimality coefficients"
187 )
188 theta = mi_master.sync_db.add_variable(
189 "theta", 0, VarType.Free, "Future profit function variable"
190 )
191 theta_fix = mi_master.sync_db.add_parameter("thetaFix", 0)
192 mi_master.instantiate(
193 "masterproblem max zmaster using lp",
194 [
195 GamsModifier(cutconst),
196 GamsModifier(cutcoeff),
197 GamsModifier(theta, UpdateAction.Fixed, theta_fix),
198 ],
199 opt,
200 )
201
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")
207 mi_sub.instantiate(
208 "subproblem max zsub using lp",
209 [GamsModifier(received), GamsModifier(demand)],
210 opt,
211 )
212
213 lower_bound = float("-inf")
214 upper_bound = float("inf")
215 obj_master = float("inf")
216 it = 1
217 theta_fix.add_record().value = 0 # fix theta for first iteration
218
219 while True:
220 print(f"Iteration: {it}")
221
222 # solve master
223 if it > 1:
224 theta_fix.clear()
225
226 mi_master.solve(SymbolUpdateType.BaseCase)
227 print(
228 f" Master {mi_master.model_status} : obj={mi_master.sync_db.get_variable('zmaster').first_record().level}"
229 )
230 if it > 1:
231 upper_bound = mi_master.sync_db.get_variable("zmaster").first_record().level
232 obj_master = (
233 mi_master.sync_db["zmaster"].first_record().level
234 - theta.first_record().level
235 )
236
237 # set received from master
238 received.clear()
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)))
242
243 cutconst.add_record(str(it))
244 obj_sub = 0.0
245 for s in data.out_db["s"]:
246 demand.clear()
247 for j in data.out_db["j"]:
248 demand.add_record(j.keys).value = scenario_data.find_record(
249 (s.key(0), j.key(0))
250 ).value
251 mi_sub.solve(SymbolUpdateType.BaseCase)
252 print(
253 f" Sub {mi_sub.model_status} : obj={mi_sub.sync_db['zsub'].first_record().level}"
254 )
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 += (
259 probability
260 * mi_sub.sync_db["market"].find_record(j.keys).marginal
261 * demand.find_record(j.keys).value
262 )
263 cutcoeff.find_record((str(it), j.key(0))).value += (
264 probability * mi_sub.sync_db["selling"].find_record(j.keys).marginal
265 )
266
267 lower_bound = max(lower_bound, obj_master + obj_sub)
268 it += 1
269 if it > max_iter:
270 raise Exception("Benders out of iterations")
271
272 print(
273 f" lower bound: {lower_bound} upper bound: {upper_bound} obj master: {obj_master}"
274 )
275 if upper_bound - lower_bound < 0.001 * (1 + abs(upper_bound)):
276 break
benders_2stage.received
GamsWorkspace received
Definition: benders_2stage.py:205
benders_2stage.cutcoeff
GamsWorkspace cutcoeff
Definition: benders_2stage.py:185
benders_2stage.cutconst
GamsWorkspace cutconst
Definition: benders_2stage.py:182
benders_2stage.demand
GamsWorkspace demand
Definition: benders_2stage.py:206