Benders2StageMT.java
1package com.gams.examples.benders;
2 
3 import java.io.File;
4 import java.util.HashMap;
5 import java.util.LinkedList;
6 import java.util.List;
7 import java.util.Map;
8 import java.util.Map.Entry;
9 
10 import com.gams.api.GAMSCheckpoint;
11 import com.gams.api.GAMSException;
12 import com.gams.api.GAMSGlobals;
13 import com.gams.api.GAMSJob;
15 import com.gams.api.GAMSModifier;
16 import com.gams.api.GAMSOptions;
17 import com.gams.api.GAMSParameter;
18 import com.gams.api.GAMSSetRecord;
19 import com.gams.api.GAMSVariable;
21 import com.gams.api.GAMSWorkspace;
23 
39 @SuppressWarnings("unchecked")
40 public class Benders2StageMT {
41 
42  public static void main(String[] args) {
43  GAMSWorkspaceInfo wsInfo = new GAMSWorkspaceInfo();
44  if (args.length > 0)
45  wsInfo.setSystemDirectory( args[0] );
46 
47  File workingDirectory = new File(System.getProperty("user.dir"), "Benders2StageMT");
48  workingDirectory.mkdir();
49  wsInfo.setWorkingDirectory(workingDirectory.getAbsolutePath());
50 
51  GAMSWorkspace ws = new GAMSWorkspace( wsInfo );
52 
53  GAMSJob dataJob = ws.addJobFromString(data);
54 
55  GAMSOptions optData = ws.addOptions();
56  optData.defines("useBig", "1");
57  optData.defines("nrScen", "100");
58 
59  dataJob.run(optData);
60 
61  optData.dispose();
62 
63  GAMSParameter scenarioData = dataJob.OutDB().getParameter("ScenarioData");
64 
65  GAMSOptions opt = ws.addOptions();
66  opt.defines("datain", dataJob.OutDB().getName());
67  int maxiter = 40;
68  opt.defines("maxiter", Integer.toString(maxiter));
69  opt.setAllModelTypes( "cplex" );
70 
71  GAMSCheckpoint cpMaster = ws.addCheckpoint();
72  GAMSCheckpoint cpSub = ws.addCheckpoint();
73 
74  ws.addJobFromString(masterModel).run(opt, cpMaster, dataJob.OutDB());
75 
76  GAMSModelInstance masteri = cpMaster.addModelInstance();
77  GAMSParameter cutconst = masteri.SyncDB().addParameter("cutconst", 1, "Benders optimality cut constant");
78  GAMSParameter cutcoeff = masteri.SyncDB().addParameter("cutcoeff", 2, "Benders optimality coefficients");
79  GAMSVariable theta = masteri.SyncDB().addVariable("theta", 0, GAMSGlobals.VarType.FREE, "Future profit function variable");
80  GAMSParameter thetaFix = masteri.SyncDB().addParameter("thetaFix", 0, "");
81  masteri.instantiate("masterproblem max zmaster using lp", opt,
82  new GAMSModifier[] { new GAMSModifier(cutconst), new GAMSModifier(cutcoeff), new GAMSModifier(theta,GAMSGlobals.UpdateAction.FIXED,thetaFix) }
83  );
84 
85  ws.addJobFromString(subModel).run(opt, cpSub, dataJob.OutDB());
86 
87 
88  int numThreads = 4;
89  GAMSModelInstance[] subi = new GAMSModelInstance[numThreads];
90  LinkedList<Tuple> demQueue = new LinkedList<Tuple>();
91 
92  subi[0] = cpSub.addModelInstance();
93  GAMSParameter received = subi[0].SyncDB().addParameter("received", 1, "units received from first stage solution");
94  GAMSParameter demand = subi[0].SyncDB().addParameter("demand", 1, "stochastic demand");
95  subi[0].instantiate("subproblem max zsub using lp", opt, new GAMSModifier[] { new GAMSModifier(received), new GAMSModifier(demand) } );
96 
97  for (int i = 1; i < numThreads; i++)
98  {
99  subi[i] = subi[0].copyModelInstance();
100  }
101 
102  opt.dispose();
103 
104  double lowerbound = Double.NEGATIVE_INFINITY, upperbound = Double.POSITIVE_INFINITY, objmaster = Double.POSITIVE_INFINITY;
105  int iter = 1;
106  do {
107  System.out.println("Iteration: " + iter);
108  // Solve master
109  if (iter == 1) // fix theta for first iteration
110  thetaFix.addRecord().setValue( 0.0 );
111  else
112  thetaFix.clear();
113 
115  System.out.println(" Master " + masteri.getModelStatus() + " : obj=" + masteri.SyncDB().getVariable("zmaster").getFirstRecord().getLevel());
116 
117  if (iter > 1)
118  upperbound = masteri.SyncDB().getVariable("zmaster").getFirstRecord().getLevel();
119 
120  objmaster = masteri.SyncDB().getVariable("zmaster").getFirstRecord().getLevel() - theta.getFirstRecord().getLevel();
121 
122  for (GAMSSetRecord s : dataJob.OutDB().getSet("s"))
123  {
124  Map<String,Double> demDict = new HashMap<String, Double>();
125  for (GAMSSetRecord j : dataJob.OutDB().getSet("j")) {
126  String[] keys = new String[] { s.getKey(0), j.getKey(0) };
127  demDict.put( j.getKey(0), new Double( scenarioData.findRecord( keys ).getValue()) );
128  }
129  String item1 = s.getKey(0);
130  Double item2 = new Double( scenarioData.findRecord(new String[] { s.getKey(0), "prob" } ).getValue() );
131  Tuple items = new Tuple( item1, item2, demDict );
132  demQueue.add(items);
133  }
134 
135  for (int i = 0; i < numThreads; i++)
136  subi[i].SyncDB().getParameter("received").clear();
137 
138  for (GAMSVariableRecord r : masteri.SyncDB().getVariable("received"))
139  {
140  cutcoeff.addRecord(new String[] {Integer.toString(iter), r.getKey(0)});
141  for (int i = 0; i < numThreads; i++)
142  {
143  subi[i].SyncDB().getParameter("received").addRecord(r.getKeys()).setValue( r.getLevel() );
144  }
145  }
146 
147  cutconst.addRecord(Integer.toString(iter));
148  double objsubsum = 0.0;
149 
150  // solve multiple model instances in parallel
151  Object queueMutex = new Object();
152  Object ioMutex = new Object();
153  Wrapper<Double>[] objsub = new Wrapper[numThreads];
154  Object[] coef = new Object[numThreads];
155  Wrapper<Double>[] cons = new Wrapper[numThreads] ;
156 
157  for (int i = 0; i < numThreads; i++)
158  {
159  objsub[i] = new Wrapper<Double>(new Double(0.0));
160  cons[i] = new Wrapper<Double>(new Double(0.0));
161  coef[i] = new HashMap<String, Double>();
162  for (GAMSSetRecord j : dataJob.OutDB().getSet("j")) {
163  Map<String, Double> cmap = (Map<String, Double>) coef[i];
164  cmap.put( j.getKey(0), new Double(0.0) );
165  }
166  }
167 
168  // Solve of sub problems in parallel
169  Scenario[] sc = new Scenario[numThreads];
170  for (int i = 0; i < numThreads; i++)
171  {
172  sc[i] = new Scenario( i, subi[i], cons[i], (Map<String, Double>) coef[i], demQueue, objsub[i], queueMutex, ioMutex );
173  sc[i].start();
174  }
175 
176 
177  // Synchronize all sub problems
178  for (int i = 0; i < numThreads; i++)
179  {
180  try {
181  sc[i].join();
182  } catch (InterruptedException e) {
183  e.printStackTrace();
184  }
185  }
186 
187  for (int i = 0; i < numThreads; i++)
188  {
189  objsubsum += objsub[i].get().doubleValue();
190  double new_consValue = cutconst.findRecord( Integer.toString(iter) ).getValue() + cons[i].get().doubleValue();
191  cutconst.findRecord( Integer.toString(iter) ).setValue( new_consValue );
192 
193  for (GAMSSetRecord j : dataJob.OutDB().getSet("j"))
194  {
195  Map<String, Double> map = (Map<String, Double>) coef[i];
196  String[] keys = new String[] { Integer.toString(iter), j.getKey(0) };
197  double newvalue = cutcoeff.findRecord( keys ).getValue( ) + map.get( j.getKey(0) ).doubleValue();
198  cutcoeff.findRecord( keys ).setValue( newvalue );
199  }
200  }
201 
202  lowerbound = Math.max(lowerbound, objmaster + objsubsum);
203 
204  iter++;
205  if (iter == maxiter + 1)
206  throw new GAMSException("Benders out of iterations");
207 
208  System.out.println(" lowerbound: " + lowerbound + " upperbound: " + upperbound + " objmaster: " + objmaster);
209 
210  } while ((upperbound - lowerbound) >= 0.001 * (1 + Math.abs(upperbound)));
211 
212  masteri.dispose();
213 
214  for(GAMSModelInstance inst : subi)
215  {
216  inst.dispose();
217  }
218  }
219 
221  static class Scenario extends Thread {
222  int _i;
223  GAMSModelInstance _submi;
224  Wrapper<Double> _cutconst;
225  List<Tuple> _demQueue;
226  Map<String, Double> _cutcoeff;
227  Wrapper<Double> _objsub;
228  Object _queueMutex;
229  Object _ioMutex;
230 
241  public Scenario(int i , GAMSModelInstance submi, Wrapper<Double> cutconst, Map<String, Double> cutcoeff, LinkedList<Tuple> demQueue,
242  Wrapper<Double> objsub, Object queueMutex, Object ioMutex) {
243  _i = i;
244  _submi = submi;
245  _cutconst = cutconst;
246  _cutcoeff = cutcoeff;
247  _demQueue = demQueue;
248  _objsub = objsub;
249  _queueMutex = queueMutex;
250  _ioMutex = ioMutex;
251  }
252 
254  public void run() {
255  while (true)
256  {
257  Tuple demDict;
258 
259  synchronized (_queueMutex)
260  {
261  if (_demQueue.size() == 0)
262  break;
263  else
264  demDict = _demQueue.remove(0); // dequeue
265  }
266 
267  _submi.SyncDB().getParameter("demand").clear();
268  for (Entry<String, Double> kv : demDict.getItem3().entrySet())
269  _submi.SyncDB().getParameter("demand").addRecord(kv.getKey()).setValue( kv.getValue() );
270 
272 
273 
274  synchronized (_ioMutex)
275  {
276  System.out.println(" Thread "+_i+" : Sub " + _submi.getModelStatus().toString() + " : obj=" + _submi.SyncDB().getVariable("zsub").getFirstRecord().getLevel());
277  }
278 
279  double probability = demDict.getItem2().doubleValue();
280  double new_objsubValue = _objsub.get().doubleValue() + probability * _submi.SyncDB().getVariable("zsub").getFirstRecord().getLevel();
281  _objsub.set( new Double( new_objsubValue ) );
282 
283  for (Entry<String, Double> kv : demDict.getItem3().entrySet())
284  {
285  double new_custconstValue = _cutconst.get().doubleValue() + probability * _submi.SyncDB().getEquation("market").findRecord( kv.getKey() ).getMarginal() * kv.getValue().doubleValue();
286  _cutconst.set( new Double( new_custconstValue ) );
287  double new_cutcoeffValue = _cutcoeff.get( kv.getKey() ).doubleValue() + probability * _submi.SyncDB().getEquation("selling").findRecord(kv.getKey()).getMarginal();
288  _cutcoeff.put( kv.getKey(), new_cutcoeffValue );
289  }
290  }
291  }
292  }
293 
294  static String data = "Sets \n"+
295  "i factories /f1*f3/ \n"+
296  "j distribution centers /d1*d5/ \n"+
297  " \n"+
298  "Parameter \n"+
299  "capacity(i) unit capacity at factories \n"+
300  " /f1 500, f2 450, f3 650/ \n"+
301  "demand(j) unit demand at distribution centers \n"+
302  " /d1 160, d2 120, d3 270, d4 325, d5 700 / \n"+
303  "prodcost unit production cost /14/ \n"+
304  "price sales price /24/ \n"+
305  "wastecost cost of removal of overstocked products /4/ \n"+
306  " \n"+
307  "Table transcost(i,j) unit transportation cost \n"+
308  " d1 d2 d3 d4 d5 \n"+
309  " f1 2.49 5.21 3.76 4.85 2.07 \n"+
310  " f2 1.46 2.54 1.83 1.86 4.76 \n"+
311  " f3 3.26 3.08 2.60 3.76 4.45; \n"+
312  " \n"+
313  "$ifthen not set useBig \n"+
314  "Set \n"+
315  " s scenarios /lo,mid,hi/ \n"+
316  " \n"+
317  "Table ScenarioData(s,*) possible outcomes for demand plus probabilities \n"+
318  " d1 d2 d3 d4 d5 prob \n"+
319  " lo 150 100 250 300 600 0.25 \n"+
320  " mid 160 120 270 325 700 0.50 \n"+
321  " hi 170 135 300 350 800 0.25; \n"+
322  "$else \n"+
323  "$if not set nrScen $set nrScen 10 \n"+
324  "Set s scenarios /s1*s%nrScen%/;\n"+
325  "parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;\n"+
326  "option seed=1234; \n"+
327  "ScenarioData(s,'prob') = 1/card(s); \n"+
328  "ScenarioData(s,j) = demand(j)*uniform(0.6,1.4); \n"+
329  "$endif \n"+
330  " \n";
331 
332  static String masterModel = "Sets \n"+
333  "i factories \n"+
334  "j distribution centers \n"+
335  " \n"+
336  "Parameter \n"+
337  "capacity(i) unit capacity at factories \n"+
338  "prodcost unit production cost \n"+
339  "transcost(i,j) unit transportation cost \n"+
340  " \n"+
341  "$if not set datain $abort 'datain not set' \n"+
342  "$gdxin %datain% \n"+
343  "$load i j capacity prodcost transcost \n"+
344  " \n"+
345  "* Benders master problem \n"+
346  "$if not set maxiter $set maxiter 25 \n"+
347  "Set \n"+
348  " iter max Benders iterations /1*%maxiter%/ \n"+
349  " \n"+
350  "Parameter \n"+
351  " cutconst(iter) constants in optimality cuts \n"+
352  " cutcoeff(iter,j) coefficients in optimality cuts \n"+
353  " \n"+
354  "Variables \n"+
355  " ship(i,j) shipments \n"+
356  " product(i) production \n"+
357  " received(j) quantity sent to market \n"+
358  " zmaster objective variable of master problem \n"+
359  " theta future profit \n"+
360  "Positive Variables ship; \n"+
361  " \n"+
362  "Equations \n"+
363  " masterobj master objective function \n"+
364  " production(i) calculate production in each factory \n"+
365  " receive(j) calculate quantity to be send to markets \n"+
366  " optcut(iter) Benders optimality cuts; \n"+
367  " \n"+
368  "masterobj.. \n"+
369  " zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j)) \n"+
370  " - sum(i,prodcost*product(i)); \n"+
371  " \n"+
372  "receive(j).. received(j) =e= sum(i, ship(i,j)); \n"+
373  " \n"+
374  "production(i).. product(i) =e= sum(j, ship(i,j)); \n"+
375  "product.up(i) = capacity(i); \n"+
376  " \n"+
377  "optcut(iter).. theta =l= cutconst(iter) + \n"+
378  " sum(j, cutcoeff(iter,j)*received(j)); \n"+
379  " \n"+
380  "model masterproblem /all/; \n"+
381  " \n"+
382  "* Initialize cut to be non-binding \n"+
383  "cutconst(iter) = 1e15; \n"+
384  "cutcoeff(iter,j) = eps; \n"+
385  " \n";
386 
387  static String subModel = "Sets \n"+
388  " i factories \n"+
389  " j distribution centers \n"+
390  " \n"+
391  "Parameter \n"+
392  " demand(j) unit demand at distribution centers \n"+
393  " price sales price \n"+
394  " wastecost cost of removal of overstocked products \n"+
395  " received(j) first stage decision units received \n"+
396  " \n"+
397  "$if not set datain $abort 'datain not set' \n"+
398  "$gdxin %datain% \n"+
399  "$load i j demand price wastecost \n"+
400  " \n"+
401  "* Benders' subproblem \n"+
402  " \n"+
403  "Variables \n"+
404  " sales(j) sales (actually sold) \n"+
405  " waste(j) overstocked products \n"+
406  " zsub objective variable of sub problem \n"+
407  "Positive variables sales, waste \n"+
408  " \n"+
409  "Equations \n"+
410  " subobj subproblem objective function \n"+
411  " selling(j) part of received is sold \n"+
412  " market(j) upperbound on sales \n"+
413  "; \n"+
414  " \n"+
415  "subobj.. \n"+
416  " zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j)); \n"+
417  " \n"+
418  "selling(j).. sales(j) + waste(j) =e= received(j); \n"+
419  " \n"+
420  "market(j).. sales(j) =l= demand(j); \n"+
421  " \n"+
422  "model subproblem /subobj,selling,market/; \n"+
423  " \n"+
424  "* Initialize received \n"+
425  "received(j) = demand(j); \n"+
426  " \n";
427 }
428 
430 class Wrapper<T> {
431  T _value;
432  public Wrapper(T value) { _value = value; }
433  public T get() { return _value; }
434  public void set(T anotherValue) { _value = anotherValue; }
435 }
436 
440 class Tuple {
441  String _item1;
442  Double _item2;
443  Map<String, Double> _item3;
444  public Tuple(String item1, Double item2, Map<String, Double> item3) {
445  _item1 = item1;
446  _item2 = item2;
447  _item3 = item3;
448  }
449  public String getItem1() { return _item1; }
450  public Double getItem2() { return _item2; }
451  public Map<String, Double> getItem3() { return _item3; }
452 }
GAMSEquation getEquation(String identifier)
GAMSGlobals.ModelStat getModelStatus()
GAMSParameter getParameter(String identifier)
void defines(String defStr, String asStr)
void setSystemDirectory(String directory)
GAMSModelInstance copyModelInstance()
GAMSParameter addParameter(String identifier, int dimension)
GAMSVariable addVariable(String identifier, int dimension, GAMSGlobals.VarType varType)
GAMSSet getSet(String identifier)
GAMSJob addJobFromString(String source)
void instantiate(String modelDefinition, GAMSModifier ... modifiers)
GAMSCheckpoint addCheckpoint()
GAMSModelInstance addModelInstance()
void setWorkingDirectory(String directory)
GAMSDatabase OutDB()
GAMSVariable getVariable(String identifier)
T addRecord(Vector< String > keys)
This example demonstrates a parallel implementation of a simple Benders decomposition method for a st...
void setAllModelTypes(String value)