Benders2Stage.cs
1using System;
2 using System.Collections.Generic;
3 using System.Text;
4 using System.IO;
5 using GAMS;
6 
8 {
25  {
26  static void Main(string[] args)
27  {
28  GAMSWorkspace ws;
29  if (Environment.GetCommandLineArgs().Length > 1)
30  ws = new GAMSWorkspace(systemDirectory: Environment.GetCommandLineArgs()[1]);
31  else
32  ws = new GAMSWorkspace();
33  GAMSJob data = ws.AddJobFromString(GetDataText());
34 
35  GAMSOptions optData = ws.AddOptions();
36  optData.Defines.Add("useBig", "1");
37  optData.Defines.Add("nrScen", "100");
38 
39  data.Run(optData);
40 
41  optData.Dispose();
42  GAMSParameter scenarioData = data.OutDB.GetParameter("ScenarioData");
43 
44  GAMSOptions opt = ws.AddOptions();
45  opt.Defines.Add("datain", data.OutDB.Name);
46  int maxiter = 40;
47  opt.Defines.Add("maxiter", maxiter.ToString());
48  opt.AllModelTypes = "cplex";
49 
50  GAMSCheckpoint cpMaster = ws.AddCheckpoint();
51  GAMSCheckpoint cpSub = ws.AddCheckpoint();
52 
53  ws.AddJobFromString(GetMasterText()).Run(opt, cpMaster, data.OutDB);
54 
55  GAMSModelInstance masteri = cpMaster.AddModelInstance();
56  GAMSParameter cutconst = masteri.SyncDB.AddParameter("cutconst", 1, "Benders optimality cut constant");
57  GAMSParameter cutcoeff = masteri.SyncDB.AddParameter("cutcoeff", 2, "Benders optimality coefficients");
58  GAMSVariable theta = masteri.SyncDB.AddVariable("theta", 0, VarType.Free, "Future profit function variable");
59  GAMSParameter thetaFix = masteri.SyncDB.AddParameter("thetaFix", 0, "");
60  masteri.Instantiate("masterproblem max zmaster using lp", opt, new GAMSModifier(cutconst), new GAMSModifier(cutcoeff), new GAMSModifier(theta, UpdateAction.Fixed,thetaFix));
61 
62  ws.AddJobFromString(GetSubText()).Run(opt, cpSub, data.OutDB);
63 
64  GAMSModelInstance subi = cpSub.AddModelInstance();
65  GAMSParameter received = subi.SyncDB.AddParameter("received", 1, "units received from master");
66  GAMSParameter demand = subi.SyncDB.AddParameter("demand", 1, "stochastic demand");
67  subi.Instantiate("subproblem max zsub using lp", opt, new GAMSModifier(received), new GAMSModifier(demand));
68 
69  opt.Dispose();
70 
71  double lowerbound = double.NegativeInfinity, upperbound = double.PositiveInfinity, objmaster = double.PositiveInfinity;
72  int iter = 1;
73  do
74  {
75  Console.WriteLine("Iteration: " + iter);
76  // Solve master
77  if (1 == iter) // fix theta for first iteration
78  thetaFix.AddRecord().Value = 0;
79  else
80  thetaFix.Clear();
81 
82  masteri.Solve(GAMSModelInstance.SymbolUpdateType.BaseCase);
83  Console.WriteLine(" Master " + masteri.ModelStatus + " : obj=" + masteri.SyncDB.GetVariable("zmaster").FirstRecord().Level);
84  if (1 < iter)
85  upperbound = masteri.SyncDB.GetVariable("zmaster").FirstRecord().Level;
86  objmaster = masteri.SyncDB.GetVariable("zmaster").FirstRecord().Level - theta.FirstRecord().Level;
87 
88  // Set received from master
89  received.Clear();
90  foreach (GAMSVariableRecord r in masteri.SyncDB.GetVariable("received"))
91  {
92  received.AddRecord(r.Keys).Value = r.Level;
93  cutcoeff.AddRecord(iter.ToString(), r.Key(0));
94  }
95 
96  cutconst.AddRecord(iter.ToString());
97  double objsub = 0.0;
98  // Benders Cut: theta <= sum(s, pi(s)*v^T(s)(b-A1*x)) v is dual soution of subproblem, A1 is first stage constraint matrix
99  foreach (GAMSSetRecord s in data.OutDB.GetSet("s"))
100  {
101  demand.Clear();
102  foreach (GAMSSetRecord j in data.OutDB.GetSet("j"))
103  demand.AddRecord(j.Keys).Value = scenarioData.FindRecord(s.Key(0), j.Key(0)).Value;
104 
106  Console.WriteLine(" Sub " + subi.ModelStatus + " : obj=" + subi.SyncDB.GetVariable("zsub").FirstRecord().Level);
107 
108 
109  double probability = scenarioData.FindRecord(s.Key(0), "prob").Value;
110  objsub += probability * subi.SyncDB.GetVariable("zsub").FirstRecord().Level;
111  foreach (GAMSSetRecord j in data.OutDB.GetSet("j"))
112  {
113  cutconst.FindRecord(iter.ToString()).Value += probability * subi.SyncDB.GetEquation("market").FindRecord(j.Keys).Marginal * demand.FindRecord(j.Keys).Value;
114  cutcoeff.FindRecord(iter.ToString(), j.Key(0)).Value += probability * subi.SyncDB.GetEquation("selling").FindRecord(j.Keys).Marginal;
115  }
116  }
117  lowerbound = Math.Max(lowerbound, objmaster + objsub);
118  iter++;
119  if (iter == maxiter + 1)
120  throw new Exception("Benders out of iterations");
121 
122  Console.WriteLine(" lowerbound: " + lowerbound + " upperbound: " + upperbound + " objmaster: " + objmaster);
123  } while ((upperbound - lowerbound) >= 0.001 * (1 + Math.Abs(upperbound)));
124 
125  masteri.Dispose();
126  subi.Dispose();
127  }
128 
129  static String GetDataText()
130  {
131  String model = @"
132 Sets
133  i factories /f1*f3/
134  j distribution centers /d1*d5/
135 
136 Parameter
137  capacity(i) unit capacity at factories
138  /f1 500, f2 450, f3 650/
139  demand(j) unit demand at distribution centers
140  /d1 160, d2 120, d3 270, d4 325, d5 700 /
141  prodcost unit production cost /14/
142  price sales price /24/
143  wastecost cost of removal of overstocked products /4/
144 
145 Table transcost(i,j) unit transportation cost
146  d1 d2 d3 d4 d5
147  f1 2.49 5.21 3.76 4.85 2.07
148  f2 1.46 2.54 1.83 1.86 4.76
149  f3 3.26 3.08 2.60 3.76 4.45;
150 
151 $ifthen not set useBig
152 Set
153  s scenarios /lo,mid,hi/
154 
155 table ScenarioData(s,*) possible outcomes for demand plus probabilities
156  d1 d2 d3 d4 d5 prob
157 lo 150 100 250 300 600 0.25
158 mid 160 120 270 325 700 0.50
159 hi 170 135 300 350 800 0.25;
160 $else
161 $if not set nrScen $set nrScen 10
162 Set s scenarios /s1*s%nrScen%/;
163 parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;
164 option seed=1234;
165 ScenarioData(s,'prob') = 1/card(s);
166 ScenarioData(s,j) = demand(j)*uniform(0.6,1.4);
167 $endif
168 ";
169 
170  return model;
171  }
172 
173  static String GetMasterText()
174  {
175  String model = @"
176 Sets
177  i factories
178  j distribution centers
179 
180 Parameter
181  capacity(i) unit capacity at factories
182  prodcost unit production cost
183  transcost(i,j) unit transportation cost
184 
185 $if not set datain $abort 'datain not set'
186 $gdxin %datain%
187 $load i j capacity prodcost transcost
188 
189 * Benders master problem
190 $if not set maxiter $set maxiter 25
191 Set
192  iter max Benders iterations /1*%maxiter%/
193 
194 Parameter
195  cutconst(iter) constants in optimality cuts
196  cutcoeff(iter,j) coefficients in optimality cuts
197 
198 Variables
199  ship(i,j) shipments
200  product(i) production
201  received(j) quantity sent to market
202  zmaster objective variable of master problem
203  theta future profit
204 Positive Variables ship;
205 
206 Equations
207  masterobj master objective function
208  production(i) calculate production in each factory
209  receive(j) calculate quantity to be send to markets
210  optcut(iter) Benders optimality cuts;
211 
212 masterobj..
213  zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j))
214  - sum(i,prodcost*product(i));
215 
216 receive(j).. received(j) =e= sum(i, ship(i,j));
217 
218 production(i).. product(i) =e= sum(j, ship(i,j));
219 product.up(i) = capacity(i);
220 
221 optcut(iter).. theta =l= cutconst(iter) +
222  sum(j, cutcoeff(iter,j)*received(j));
223 
224 model masterproblem /all/;
225 
226 * Initialize cut to be non-binding
227 cutconst(iter) = 1e15;
228 cutcoeff(iter,j) = eps;
229 ";
230 
231  return model;
232  }
233 
234  static String GetSubText()
235  {
236  String model = @"
237 Sets
238  i factories
239  j distribution centers
240 
241 Parameter
242  demand(j) unit demand at distribution centers
243  price sales price
244  wastecost cost of removal of overstocked products
245  received(j) first stage decision units received
246 
247 $if not set datain $abort 'datain not set'
248 $gdxin %datain%
249 $load i j demand price wastecost
250 
251 * Benders' subproblem
252 
253 Variables
254  sales(j) sales (actually sold)
255  waste(j) overstocked products
256  zsub objective variable of sub problem
257 Positive variables sales, waste
258 
259 Equations
260  subobj subproblem objective function
261  selling(j) part of received is sold
262  market(j) upperbound on sales
263 ;
264 
265 subobj..
266  zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
267 
268 selling(j).. sales(j) + waste(j) =e= received(j);
269 
270 market(j).. sales(j) =l= demand(j);
271 
272 model subproblem /subobj,selling,market/;
273 
274 * Initialize received
275 received(j) = demand(j);
276 ";
277 
278  return model;
279  }
280 
281  }
282 }
This example demonstrates a sequential implementation of a simple Benders decomposition method for a ...
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)
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)