Benders2StageMT.cs
1using System;
2 using System.Collections.Generic;
3 using System.Text;
4 using System.IO;
5 using System.Threading.Tasks;
6 using 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  {
62  GAMSWorkspace ws;
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 = @"
189 Sets
190  i factories /f1*f3/
191  j distribution centers /d1*d5/
192 
193 Parameter
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 
202 Table 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
209 Set
210  s scenarios /lo,mid,hi/
211 
212 table ScenarioData(s,*) possible outcomes for demand plus probabilities
213  d1 d2 d3 d4 d5 prob
214 lo 150 100 250 300 600 0.25
215 mid 160 120 270 325 700 0.50
216 hi 170 135 300 350 800 0.25;
217 $else
218 $if not set nrScen $set nrScen 10
219 Set s scenarios /s1*s%nrScen%/;
220 parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;
221 option seed=1234;
222 ScenarioData(s,'prob') = 1/card(s);
223 ScenarioData(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 = @"
233 Sets
234  i factories
235  j distribution centers
236 
237 Parameter
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
248 Set
249  iter max Benders iterations /1*%maxiter%/
250 
251 Parameter
252  cutconst(iter) constants in optimality cuts
253  cutcoeff(iter,j) coefficients in optimality cuts
254 
255 Variables
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
261 Positive Variables ship;
262 
263 Equations
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 
269 masterobj..
270  zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j))
271  - sum(i,prodcost*product(i));
272 
273 receive(j).. received(j) =e= sum(i, ship(i,j));
274 
275 production(i).. product(i) =e= sum(j, ship(i,j));
276 product.up(i) = capacity(i);
277 
278 optcut(iter).. theta =l= cutconst(iter) +
279  sum(j, cutcoeff(iter,j)*received(j));
280 
281 model masterproblem /all/;
282 
283 * Initialize cut to be non-binding
284 cutconst(iter) = 1e15;
285 cutcoeff(iter,j) = eps;
286 ";
287 
288  return model;
289  }
290 
291  static String GetSubText()
292  {
293  String model = @"
294 Sets
295  i factories
296  j distribution centers
297 
298 Parameter
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 
310 Variables
311  sales(j) sales (actually sold)
312  waste(j) overstocked products
313  zsub objective variable of sub problem
314 Positive variables sales, waste
315 
316 Equations
317  subobj subproblem objective function
318  selling(j) part of received is sold
319  market(j) upperbound on sales
320 ;
321 
322 subobj..
323  zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
324 
325 selling(j).. sales(j) + waste(j) =e= received(j);
326 
327 market(j).. sales(j) =l= demand(j);
328 
329 model subproblem /subobj,selling,market/;
330 
331 * Initialize received
332 received(j) = demand(j);
333 ";
334 
335  return model;
336  }
337 
338  }
339 }
void Instantiate(string modelDefinition, params GAMSModifier[] modifiers)
string Key(int index)
Dictionary< string, string > Defines
new GAMSParameterRecord AddRecord(params string[] keys)
GAMSCheckpoint AddCheckpoint(string checkpointName=null)
GAMSDatabase OutDB
GAMSParameter GetParameter(string parameterIdentifier)
This example demonstrates a parallel implementation of a simple Benders decomposition method for a st...
GAMSOptions AddOptions(GAMSOptions optFrom=null)
UpdateAction
GAMSModelInstance AddModelInstance(string modelInstanceName=null)
new GAMSVariableRecord FirstRecord()
GAMSSet GetSet(string setIdentifier)
GAMSVariable GetVariable(string variableIdentifier)
new GAMSEquationRecord FindRecord(params string[] keys)
new GAMSParameterRecord FindRecord(params string[] keys)
GAMSJob AddJobFromString(string gamsSource, GAMSCheckpoint checkpoint=null, string jobName=null)
void Solve(SymbolUpdateType updateType=SymbolUpdateType.BaseCase, TextWriter output=null, GAMSModelInstanceOpt miOpt=null)
GAMSParameter AddParameter(string identifier, int dimension, string explanatoryText="")
GAMSVariable AddVariable(string identifier, int dimension, VarType varType, string explanatoryText="")
void Run(GAMSOptions gamsOptions=null, GAMSCheckpoint checkpoint=null, TextWriter output=null, Boolean createOutDB=true)
GAMSEquation GetEquation(string equationIdentifier)