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