Loading...
Searching...
No Matches
Benders2StageMT.cs
1using System;
2using System.Collections.Generic;
3using System.Text;
4using System.IO;
5using System.Threading.Tasks;
6using GAMS;
7
9{
27 {
28 private static void ScenSolve(GAMSModelInstance subi, ref double cutconst, ref Dictionary<string, double> cutcoeff, Queue<Tuple<string, double, Dictionary<string, double>>> demQueue, ref double objsub, Object queueMutex, Object ioMutex)
29 {
30 while (true)
31 {
32 Tuple<string, double, Dictionary<string, double>> demDict;
33 lock (queueMutex)
34 {
35 if (0 == demQueue.Count)
36 return;
37 demDict = demQueue.Dequeue();
38 }
39
40 subi.SyncDB.GetParameter("demand").Clear();
41
42 foreach (KeyValuePair<string, double> kv in demDict.Item3)
43 subi.SyncDB.GetParameter("demand").AddRecord(kv.Key).Value = kv.Value;
44
46
47 lock (ioMutex)
48 Console.WriteLine(" Sub " + subi.ModelStatus + " : obj=" + subi.SyncDB.GetVariable("zsub").FirstRecord().Level);
49
50 double probability = demDict.Item2;
51 objsub += probability * subi.SyncDB.GetVariable("zsub").FirstRecord().Level;
52
53 foreach (KeyValuePair<string, double> kv in demDict.Item3)
54 {
55 cutconst += probability * subi.SyncDB.GetEquation("market").FindRecord(kv.Key).Marginal * kv.Value;
56 cutcoeff[kv.Key] += probability * subi.SyncDB.GetEquation("selling").FindRecord(kv.Key).Marginal;
57 }
58 }
59 }
60
61 static void Main(string[] args)
62 {
64 if (Environment.GetCommandLineArgs().Length > 1)
65 ws = new GAMSWorkspace(systemDirectory: Environment.GetCommandLineArgs()[1]);
66 else
67 ws = new GAMSWorkspace();
68 GAMSJob data = ws.AddJobFromString(GetDataText());
69
70 GAMSOptions optData = ws.AddOptions();
71 optData.Defines.Add("useBig", "1");
72 optData.Defines.Add("nrScen", "100");
73
74 data.Run(optData);
75
76 optData.Dispose();
77 GAMSParameter scenarioData = data.OutDB.GetParameter("ScenarioData");
78
79 GAMSOptions opt = ws.AddOptions();
80 opt.Defines.Add("datain", data.OutDB.Name);
81 int maxiter = 40;
82 opt.Defines.Add("maxiter", maxiter.ToString());
83 opt.AllModelTypes = "cplex";
84
85 GAMSCheckpoint cpMaster = ws.AddCheckpoint();
86 GAMSCheckpoint cpSub = ws.AddCheckpoint();
87
88 ws.AddJobFromString(GetMasterText()).Run(opt, cpMaster, data.OutDB);
89
90 GAMSModelInstance masteri = cpMaster.AddModelInstance();
91 GAMSParameter cutconst = masteri.SyncDB.AddParameter("cutconst", 1, "Benders optimality cut constant");
92 GAMSParameter cutcoeff = masteri.SyncDB.AddParameter("cutcoeff", 2, "Benders optimality coefficients");
93 GAMSVariable theta = masteri.SyncDB.AddVariable("theta", 0, VarType.Free, "Future profit function variable");
94 GAMSParameter thetaFix = masteri.SyncDB.AddParameter("thetaFix", 0, "");
95 masteri.Instantiate("masterproblem max zmaster using lp", opt, new GAMSModifier(cutconst), new GAMSModifier(cutcoeff), new GAMSModifier(theta, UpdateAction.Fixed, thetaFix));
96
97 ws.AddJobFromString(GetSubText()).Run(opt, cpSub, data.OutDB);
98
99 int numThreads = 4;
100 GAMSModelInstance[] subi = new GAMSModelInstance[numThreads];
101 Queue<Tuple<string, double, Dictionary<string, double>>> demQueue = new Queue<Tuple<string, double, Dictionary<string, double>>>();
102
103 subi[0] = cpSub.AddModelInstance();
104 GAMSParameter received = subi[0].SyncDB.AddParameter("received", 1, "units received from first stage solution");
105 GAMSParameter demand = subi[0].SyncDB.AddParameter("demand", 1, "stochastic demand");
106 subi[0].Instantiate("subproblem max zsub using lp", opt, new GAMSModifier(received), new GAMSModifier(demand));
107
108 for (int i = 1; i < numThreads; i++)
109 subi[i] = subi[0].CopyModelinstance();
110
111 opt.Dispose();
112
113 double lowerbound = double.NegativeInfinity, upperbound = double.PositiveInfinity, objmaster = double.PositiveInfinity;
114 int iter = 1;
115 do
116 {
117 Console.WriteLine("Iteration: " + iter);
118 // Solve master
119 if (1 == iter) // fix theta for first iteration
120 thetaFix.AddRecord().Value = 0;
121 else
122 thetaFix.Clear();
123
124 masteri.Solve(GAMSModelInstance.SymbolUpdateType.BaseCase);
125 Console.WriteLine(" Master " + masteri.ModelStatus + " : obj=" + masteri.SyncDB.GetVariable("zmaster").FirstRecord().Level);
126 if (1 < iter)
127 upperbound = masteri.SyncDB.GetVariable("zmaster").FirstRecord().Level;
128 objmaster = masteri.SyncDB.GetVariable("zmaster").FirstRecord().Level - theta.FirstRecord().Level;
129
130 foreach (GAMSSetRecord s in data.OutDB.GetSet("s"))
131 {
132 Dictionary<string, double> demDict = new Dictionary<string, double>();
133 foreach (GAMSSetRecord j in data.OutDB.GetSet("j"))
134 demDict[j.Key(0)] = scenarioData.FindRecord(s.Key(0), j.Key(0)).Value;
135 demQueue.Enqueue(new Tuple<string, double, Dictionary<string, double>>(s.Key(0), scenarioData.FindRecord(s.Key(0), "prob").Value, demDict));
136 }
137
138 for (int i = 0; i < numThreads; i++)
139 subi[i].SyncDB.GetParameter("received").Clear();
140 foreach (GAMSVariableRecord r in masteri.SyncDB.GetVariable("received"))
141 {
142 cutcoeff.AddRecord(iter.ToString(), r.Key(0));
143 for (int i = 0; i < numThreads; i++)
144 subi[i].SyncDB.GetParameter("received").AddRecord(r.Keys).Value = r.Level;
145 }
146
147 cutconst.AddRecord(iter.ToString());
148 double objsubsum = 0.0;
149
150 // solve multiple model instances in parallel
151 Object queueMutex = new Object();
152 Object ioMutex = new Object();
153 double[] objsub = new double[numThreads];
154 Dictionary<string, double>[] coef = new Dictionary<string, double>[numThreads];
155 double[] cons = new double[numThreads];
156
157 for (int i = 0; i < numThreads; i++)
158 {
159 coef[i] = new Dictionary<string, double>();
160 foreach (GAMSSetRecord j in data.OutDB.GetSet("j"))
161 coef[i].Add(j.Key(0), 0.0);
162 }
163
164 Parallel.For(0, numThreads, delegate(int i) { ScenSolve(subi[i], ref cons[i], ref coef[i], demQueue, ref objsub[i], queueMutex, ioMutex); });
165
166 for (int i = 0; i < numThreads; i++)
167 {
168 objsubsum += objsub[i];
169 cutconst.FindRecord(iter.ToString()).Value += cons[i];
170 foreach (GAMSSetRecord j in data.OutDB.GetSet("j"))
171 cutcoeff.FindRecord(iter.ToString(), j.Key(0)).Value += coef[i][j.Key(0)];
172 }
173 lowerbound = Math.Max(lowerbound, objmaster + objsubsum);
174
175 iter++;
176 if (iter == maxiter + 1)
177 throw new Exception("Benders out of iterations");
178
179 Console.WriteLine(" lowerbound: " + lowerbound + " upperbound: " + upperbound + " objmaster: " + objmaster);
180 } while ((upperbound - lowerbound) >= 0.001 * (1 + Math.Abs(upperbound)));
181
182 masteri.Dispose();
183 foreach (GAMSModelInstance inst in subi)
184 inst.Dispose();
185 }
186
187 static String GetDataText()
188 {
189 String model = @"
190Sets
191 i factories /f1*f3/
192 j distribution centers /d1*d5/
193
194Parameter
195 capacity(i) unit capacity at factories
196 /f1 500, f2 450, f3 650/
197 demand(j) unit demand at distribution centers
198 /d1 160, d2 120, d3 270, d4 325, d5 700 /
199 prodcost unit production cost /14/
200 price sales price /24/
201 wastecost cost of removal of overstocked products /4/
202
203Table transcost(i,j) unit transportation cost
204 d1 d2 d3 d4 d5
205 f1 2.49 5.21 3.76 4.85 2.07
206 f2 1.46 2.54 1.83 1.86 4.76
207 f3 3.26 3.08 2.60 3.76 4.45;
208
209$ifthen not set useBig
210Set
211 s scenarios /lo,mid,hi/
212
213table ScenarioData(s,*) possible outcomes for demand plus probabilities
214 d1 d2 d3 d4 d5 prob
215lo 150 100 250 300 600 0.25
216mid 160 120 270 325 700 0.50
217hi 170 135 300 350 800 0.25;
218$else
219$if not set nrScen $set nrScen 10
220Set s scenarios /s1*s%nrScen%/;
221parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;
222option seed=1234;
223ScenarioData(s,'prob') = 1/card(s);
224ScenarioData(s,j) = demand(j)*uniform(0.6,1.4);
225$endif
226";
227
228 return model;
229 }
230
231 static String GetMasterText()
232 {
233 String model = @"
234Sets
235 i factories
236 j distribution centers
237
238Parameter
239 capacity(i) unit capacity at factories
240 prodcost unit production cost
241 transcost(i,j) unit transportation cost
242
243$if not set datain $abort 'datain not set'
244$gdxin %datain%
245$load i j capacity prodcost transcost
246
247* Benders master problem
248$if not set maxiter $set maxiter 25
249Set
250 iter max Benders iterations /1*%maxiter%/
251
252Parameter
253 cutconst(iter) constants in optimality cuts
254 cutcoeff(iter,j) coefficients in optimality cuts
255
256Variables
257 ship(i,j) shipments
258 product(i) production
259 received(j) quantity sent to market
260 zmaster objective variable of master problem
261 theta future profit
262Positive Variables ship;
263
264Equations
265 masterobj master objective function
266 production(i) calculate production in each factory
267 receive(j) calculate quantity to be send to markets
268 optcut(iter) Benders optimality cuts;
269
270masterobj..
271 zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j))
272 - sum(i,prodcost*product(i));
273
274receive(j).. received(j) =e= sum(i, ship(i,j));
275
276production(i).. product(i) =e= sum(j, ship(i,j));
277product.up(i) = capacity(i);
278
279optcut(iter).. theta =l= cutconst(iter) +
280 sum(j, cutcoeff(iter,j)*received(j));
281
282model masterproblem /all/;
283
284* Initialize cut to be non-binding
285cutconst(iter) = 1e15;
286cutcoeff(iter,j) = eps;
287";
288
289 return model;
290 }
291
292 static String GetSubText()
293 {
294 String model = @"
295Sets
296 i factories
297 j distribution centers
298
299Parameter
300 demand(j) unit demand at distribution centers
301 price sales price
302 wastecost cost of removal of overstocked products
303 received(j) first stage decision units received
304
305$if not set datain $abort 'datain not set'
306$gdxin %datain%
307$load i j demand price wastecost
308
309* Benders' subproblem
310
311Variables
312 sales(j) sales (actually sold)
313 waste(j) overstocked products
314 zsub objective variable of sub problem
315Positive variables sales, waste
316
317Equations
318 subobj subproblem objective function
319 selling(j) part of received is sold
320 market(j) upperbound on sales
321;
322
323subobj..
324 zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
325
326selling(j).. sales(j) + waste(j) =e= received(j);
327
328market(j).. sales(j) =l= demand(j);
329
330model subproblem /subobj,selling,market/;
331
332* Initialize received
333received(j) = demand(j);
334";
335
336 return model;
337 }
338
339 }
340}
This example demonstrates a parallel implementation of a simple Benders decomposition method for a st...
GAMSModelInstance AddModelInstance(string modelInstanceName=null)
GAMSVariable GetVariable(string variableIdentifier)
GAMSParameter GetParameter(string parameterIdentifier)
GAMSSet GetSet(string setIdentifier)
GAMSVariable AddVariable(string identifier, int dimension, VarType varType, string explanatoryText="")
GAMSEquation GetEquation(string equationIdentifier)
GAMSParameter AddParameter(string identifier, int dimension, string explanatoryText="")
new GAMSEquationRecord FindRecord(params string[] keys)
GAMSDatabase OutDB
void Run(GAMSOptions gamsOptions=null, GAMSCheckpoint checkpoint=null, TextWriter output=null, Boolean createOutDB=true)
void Solve(SymbolUpdateType updateType=SymbolUpdateType.BaseCase, TextWriter output=null, GAMSModelInstanceOpt miOpt=null)
void Instantiate(string modelDefinition, params GAMSModifier[] modifiers)
Dictionary< string, string > Defines
new GAMSParameterRecord FindRecord(params string[] keys)
new GAMSParameterRecord AddRecord(params string[] keys)
string Key(int index)
new GAMSVariableRecord FirstRecord()
GAMSJob AddJobFromString(string gamsSource, GAMSCheckpoint checkpoint=null, string jobName=null)
GAMSCheckpoint AddCheckpoint(string checkpointName=null)
GAMSOptions AddOptions(GAMSOptions optFrom=null)
UpdateAction