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
8namespace BendersDecomposition2StageSP
9{
26 class Benders2StageMT
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
45 subi.Solve(GAMSModelInstance.SymbolUpdateType.BaseCase);
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 {
63 GAMSWorkspace ws;
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}
BendersDecomposition2StageSP.Benders2StageMT
This example demonstrates a parallel implementation of a simple Benders decomposition method for a st...
Definition: Benders2StageMT.cs:27
GAMS::GAMSCheckpoint::AddModelInstance
GAMSModelInstance AddModelInstance(string modelInstanceName=null)
GAMS::GAMSDatabase::GetVariable
GAMSVariable GetVariable(string variableIdentifier)
GAMS::GAMSDatabase::GetParameter
GAMSParameter GetParameter(string parameterIdentifier)
GAMS::GAMSDatabase::GetSet
GAMSSet GetSet(string setIdentifier)
GAMS::GAMSDatabase::AddVariable
GAMSVariable AddVariable(string identifier, int dimension, VarType varType, string explanatoryText="")
GAMS::GAMSDatabase::GetEquation
GAMSEquation GetEquation(string equationIdentifier)
GAMS::GAMSDatabase::AddParameter
GAMSParameter AddParameter(string identifier, int dimension, string explanatoryText="")
GAMS::GAMSEquation::FindRecord
new GAMSEquationRecord FindRecord(params string[] keys)
GAMS::GAMSJob::OutDB
GAMSDatabase OutDB
GAMS::GAMSJob::Run
void Run(GAMSOptions gamsOptions=null, GAMSCheckpoint checkpoint=null, TextWriter output=null, Boolean createOutDB=true)
GAMS::GAMSModelInstance::Solve
void Solve(SymbolUpdateType updateType=SymbolUpdateType.BaseCase, TextWriter output=null, GAMSModelInstanceOpt miOpt=null)
GAMS::GAMSModelInstance::Instantiate
void Instantiate(string modelDefinition, params GAMSModifier[] modifiers)
GAMS::GAMSOptions::Defines
Dictionary< string, string > Defines
GAMS::GAMSParameter::FindRecord
new GAMSParameterRecord FindRecord(params string[] keys)
GAMS::GAMSParameter::AddRecord
new GAMSParameterRecord AddRecord(params string[] keys)
GAMS::GAMSSymbolRecord::Key
string Key(int index)
GAMS::GAMSVariable::FirstRecord
new GAMSVariableRecord FirstRecord()
GAMS::GAMSWorkspace::AddJobFromString
GAMSJob AddJobFromString(string gamsSource, GAMSCheckpoint checkpoint=null, string jobName=null)
GAMS::GAMSWorkspace::AddCheckpoint
GAMSCheckpoint AddCheckpoint(string checkpointName=null)
GAMS::GAMSWorkspace::AddOptions
GAMSOptions AddOptions(GAMSOptions optFrom=null)
GAMS::UpdateAction
UpdateAction