Loading...
Searching...
No Matches
Benders2Stage.cs
1using System;
2using System.Collections.Generic;
3using System.Text;
4using System.IO;
5using GAMS;
6
7namespace BendersDecomposition2StageSP
8{
24 class Benders2Stage
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
105 subi.Solve(GAMSModelInstance.SymbolUpdateType.BaseCase);
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 = @"
132Sets
133 i factories /f1*f3/
134 j distribution centers /d1*d5/
135
136Parameter
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
145Table 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
152Set
153 s scenarios /lo,mid,hi/
154
155table ScenarioData(s,*) possible outcomes for demand plus probabilities
156 d1 d2 d3 d4 d5 prob
157lo 150 100 250 300 600 0.25
158mid 160 120 270 325 700 0.50
159hi 170 135 300 350 800 0.25;
160$else
161$if not set nrScen $set nrScen 10
162Set s scenarios /s1*s%nrScen%/;
163parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;
164option seed=1234;
165ScenarioData(s,'prob') = 1/card(s);
166ScenarioData(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 = @"
176Sets
177 i factories
178 j distribution centers
179
180Parameter
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
191Set
192 iter max Benders iterations /1*%maxiter%/
193
194Parameter
195 cutconst(iter) constants in optimality cuts
196 cutcoeff(iter,j) coefficients in optimality cuts
197
198Variables
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
204Positive Variables ship;
205
206Equations
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
212masterobj..
213 zmaster =e= theta -sum((i,j), transcost(i,j)*ship(i,j))
214 - sum(i,prodcost*product(i));
215
216receive(j).. received(j) =e= sum(i, ship(i,j));
217
218production(i).. product(i) =e= sum(j, ship(i,j));
219product.up(i) = capacity(i);
220
221optcut(iter).. theta =l= cutconst(iter) +
222 sum(j, cutcoeff(iter,j)*received(j));
223
224model masterproblem /all/;
225
226* Initialize cut to be non-binding
227cutconst(iter) = 1e15;
228cutcoeff(iter,j) = eps;
229";
230
231 return model;
232 }
233
234 static String GetSubText()
235 {
236 String model = @"
237Sets
238 i factories
239 j distribution centers
240
241Parameter
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
253Variables
254 sales(j) sales (actually sold)
255 waste(j) overstocked products
256 zsub objective variable of sub problem
257Positive variables sales, waste
258
259Equations
260 subobj subproblem objective function
261 selling(j) part of received is sold
262 market(j) upperbound on sales
263;
264
265subobj..
266 zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));
267
268selling(j).. sales(j) + waste(j) =e= received(j);
269
270market(j).. sales(j) =l= demand(j);
271
272model subproblem /subobj,selling,market/;
273
274* Initialize received
275received(j) = demand(j);
276";
277
278 return model;
279 }
280
281 }
282}
BendersDecomposition2StageSP.Benders2Stage
This example demonstrates a sequential implementation of a simple Benders decomposition method for a ...
Definition: Benders2Stage.cs:25
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