Issue with running SDDP example available online

[priyanka_shinde]

I tried to run the SDDP example code that is available here:
https://www.gams.com/latest/gamslib_ml/ ... _sddp.html
But it says "Starting execution - empty program" and the status shows as "Normal completion". This appears in a few seconds.
Also, the lst file doesn't give any results.

On the contrary when I run the example given at http://www.gams.com/modlib/adddocs/sddp_trad.gms with traditional loops for the same problem. It takes around 2 hours to get the results but at least it works.

Can someone tell me what I might be doing wrong while trying the first code?

Also, since there is no other example GAMS code for SDDP available online, as a beginner in this topic, I will be grateful if someone is willing to share an example.

Thanks
Priyanka

[bussieck]

Most likely you did not start it with one of the LP solvers that are "required" by the sddp model:

Code: Select all

$ifI %gams.lp%==cplexd $goto continue
$ifI %gams.lp%==xpress $goto continue
$ifI %gams.lp%==gurobi $goto continue
$log SDDP only tested for selected LP solvers
$exit

You can just remove the code and then it will work with any LP solver (that is capable of solvelink=5).

-Michael

[priyanka_shinde]

Thanks for your reply Michael.
I tried to do as you suggested. Then I tried to run the rest of the code and I got these two errors:
36 '=' or '..' or ':=' or '$=' operator expected
rest of statement ignored
140 Unknown symbol

I am not sure what is wrong? Was there still a problem with choosing an LP solver?
I can post the rest of the code if that is required. Let me know.

-Priyanka

[bussieck]

Hard to say. Upload the entire lst file. Perhaps you are using an incompatible, i.e. old, GAMS version.

-Michael

[priyanka_shinde]

Hi Micheal,
I have attached the log file and here is the content from lst file:

GAMS 24.7.3 r58181 Released Jul 11, 2016 WEX-WEI x86 64bit/MS Windows 04/22/20 09:03:31 Page 1
Multi-Stage Stochastic Water Reservoir Model solved with SDDP (SDDP,SEQ=357)
C o m p i l a t i o n


2
The Stochastic Dual Dynamic Programming (SDDP) algorithm for solving
multi-stochastic linear programs uses, similar to the well known
Benders decomposition, the concept of a future cost function (FCF).
The algorithm works with an underestimate approximation of this FCF by
iterativley adding supporting hyperplanes (Benders cuts) and
therefore improving the approximation.

This implementation uses Gather-Update-Solve-Scatter (GUSS) in GAMS
making the many related solves efficient. An implementation using
traditional loops can be found at:

http://www.gams.com/modlib/adddocs/sddp_trad.gms

This SDDP algorithm has been implemented using IBM Ilog's Concert
Technology (Version 12.2). The C++ source code is available at

http://www.gams.com/modlib/adddocs/sddp.cpp

The data and core model has been provided by Vattenfall Energy Trading.


Pereira, M V F, and Pinto, L M V G, Stochastic optimization of a
multireservoir hydroelectric system: A decomposition approach. Water
Resources Research 21, 6 (1985), 779-792.

Pereira, M V F, and Pinto, L M V G, Multi-stage stochastic
optimization applied to energy planning. Mathematical Programming 52
(1991), 359-375.

Keywords: linear programming, stochastic programming, scenario analysis, G
USS,
water resource management
36
38
39 * Only test fast LP
40
41
43
44 Set
45 w 'weeks' / w1*w52 /
46 t 'hours of the year' / t1*t8736 /
47 ft 'types of fuel available (plants)' / Hydro, HardCoal, Nuclear /;
48
49 Parameter
50 demand(t) 'power demand MW'
51 exchange(t) 'cross border flow MW'
52 wCapacity(w,ft) 'capacity of plant MW'
53 wPrices(w,*) 'coal and CO2 price EUR per ton'
54 wInflow(w) 'inflows into the reservoir MWh per week'
55 resmax(t) 'max reservoir level MW' / t1*t8735 106.2e6, t8736 87e6
/
56 resmin(t) 'min reservoir level MW' / t1*t8735 10e6, t8736 87e6
/;
57
66
71
GDXIN \\ug.kth.se\dfs\home\p\v\pvshinde\appdata\xp.V2\Documents\gamsdir\projdi
r\water.gdx
--- LOAD demand = 1:demand
--- LOAD exchange = 2:exchange
--- LOAD wCapacity = 3:wCapacity
--- LOAD wPrices = 4:wPrices
--- LOAD wInflow = 5:wInflow
74
75 * Week hour mapping
76 Set
77 h 'hours of one week' / h1*h168 /
78 wt(w,t) 'map weeks and hours'
79 last(w,t) 'last hour t of each week w'
80 first(w,t) 'first hour t of each week w';
81
82 loop(w,
83 loop((h,t)$sameas('t1',t), wt(w,t + ((ord(w) - 1)*card(h) + ord(h) - 1)
) = yes;);
84 loop(t$sameas('t1',t),
85 first(w,t+((ord(w) - 1)*card(h))) = yes;
86 last(w,t+(ord(w)*card(h) - 1)) = yes;
87 );
88 );
89
90 Parameter
91 capacity(t,ft) 'capacity of plant MW'
92 gencost(t,ft) 'generating cost EUR per MW' / #t.Hydro 0, #t.Nuclear 15
/
93 inflow(t) 'inflows into the reservoir MW';
94
95 inflow(t) = sum(wt(w,t), wInflow(w)/card(h));
96 capacity(t,ft) = sum(wt(w,t), wCapacity(w,ft));
97 gencost(t,'HardCoal') = sum(wt(w,t), (wPrices(w,'HardCoal') + 2.361*wPrice
s(w,'CO2')))/(6.98*0.39);
98
99 Set
100 tt(t) 'control set for hours'
101 ww(w) 'control set for weeks';
102
103 Positive Variable
104 GAP(t) 'gap generation in hour t MW'
105 RES(t) 'reservoir energy (level) at end of t MW'
106 SPILL(t) 'spillage in hour t MW'
107 X(t,ft) 'generated power in hour t by plant with fuel type MW'
108 SLACKUP(t) 'slack variable for the upper reservoir level equation MW'
109 SLACKLO(t) 'slack variable for the lower reservoir level equation MW';
110
111 Variable COST 'total generation cost EUR';
112
113 Equation
114 Cont(t) 'hydraulic continuity equation'
115 DemSat(t) 'demand satisfaction'
116 ResUp(t) 'maximum reservoir level'
117 ResLo(t) 'minimum reservoir level'
118 PlantCap(t,ft) 'plant power generation capacity'
119 Obj 'objective function';
120
121 Cont(tt(t))..
122 RES(t) - RES(t--1) + X(t,'Hydro') + SPILL(t) =e= inflow(t);
123
124 DemSat(tt(t))..
125 sum(ft$capacity(t,ft), X(t,ft)) + GAP(t) =g= demand(t) - exchange(t);
126
127 ResUp(tt(t))..
128 -RES(t) + SLACKUP(t) =g= -resmax(t);
129
130 ResLo(tt(t))..
131 RES(t) + SLACKLO(t) =g= resmin(t);
132
133 PlantCap(tt(t),ft)$capacity(t,ft)..
134 -X(t,ft) =g= -capacity(t,ft);
135
136 Obj.. COST =e= sum((tt(t),ft)$capacity(t,ft), gencost(t,ft)*X(t,ft))
137 + sum(tt(t), 1000*GAP(t) + 10e6*(SLACKUP(t) + SLACKLO(t)));
138
139 Model water / all /;
140 * Deterministic model
141 tt(t) = yes;
143
GAMS 24.7.3 r58181 Released Jul 11, 2016 WEX-WEI x86 64bit/MS Windows 04/22/20 09:03:31 Page 2
SDDP Algorithm
C o m p i l a t i o n


145
146 * Scenario data
147 Set s 'scenarios' / s1*s12 /;
148
149 Parameter sw_Inflow(w,s) 'inflow realizations MWh per week';
--- LOAD sw_Inflow = 6:sw_Inflow
151
154
155 Set
156 i 'trial solutions' / i1*i5 /
157 j 'iteration index' / j1*j20 /
158 jj(j) 'dynamic j';
159
160 Parameter
161 cont_m(j,i,w) 'dual variables associated with the mass balance constrai
nt'
162 delta(j,i,w) 'RHS of the Benders cuts';
163
164 cont_m(j,i,w) = 0;
165 delta(j,i,w) = 0;
166 jj(j) = no;
167
168 Free Variable ACOST 'approximation of cost';
169 Positive Variable ALPHA(t) 'approximation of future cost function (FCF)';
170
171 Equation
172 Obj_Approx 'objective function for the one-stage dispatch model'
173 Cuts(j,i,w,t) 'renders cuts';
174
175 Obj_Approx..
176 ACOST =e= sum((tt(t),ft)$capacity(t,ft), gencost(t,ft)*X(t,ft))
177 + sum(tt(t), 1000*GAP(t) + 10e6*(SLACKUP(t) + SLACKLO(t)))
178 + sum(last(ww,tt(t)), ALPHA(t+1));
179
180 Cuts(jj,i,last(ww(w),t))$(ord(w) < card(w))..
181 ALPHA(t+1) - cont_m(jj,i,w+1)*RES(t) =g= delta(jj,i,w);
182
183 Model waterSDDP / water, obj_approx, cuts /;
184 * No superflous output and in core solving
185 option limRow = 0, limCol = 0, solPrint = silent, solveLink = 5;
186
187 Parameter i_res(i,t) 'trial decisions';
188
189 Set jwis(j,w,i,s) 'scenario trial sample';
190
191 * Select some for the first backward recursion;
192 loop(i, i_res(i,t) = resmin(t) + (ord(i) - 1)/(card(i) - 1)*(resmax(t) - r
esmin(t)));
193 loop(s$sameas('s1',s), jwis(j,w,i,s + UniFormInt(0,card(s) - 1))$(not same
as('w1',w)) = yes);
194
195 * We need to iterate w from back to front
196 Set revt(w,w) 'backward loop for backward recursion excluding week 1';
197 loop(w$(ord(w) < card(w)), revt(w,w + (card(w) - 2*ord(w) + 1)) = yes);
198
200 Parameter
201 conv(j,*) 'convergence parameters'
202 zt(j,i) 'objective for trial solution';
203
205
220
221 Alias (w,wp,wloop), (i,ip);
222
223 Set is(i,s) / #i.#s /;
224
225 Parameter
226 is_RES(i,s,t)
227 is_inflow(i,s,t)
228 is_Cont(i,s,t)
229 is_PlantCap(i,s,t,ft)
230 is_ResUp(i,s,t)
231 is_ResLo(i,s,t)
232 is_DemSat(i,s,t)
233 is_Cuts(i,s,j,i,w,t)
234 is_ACOST(i,s)
235 is_ALPHA(i,s,t)
237 so / SkipBaseCase 1, LogOption 1, UpdateType 2 /;
238
239 Set
240 dict_b 'backward recursion'
241 / is .scenario.''
242 so .opt .''
243 RES .fixed .is_RES
244 inflow .param .is_inflow
245 Cont .marginal.is_Cont
246 PlantCap.marginal.is_PlantCap
247 ResUp .marginal.is_ResUp
248 ResLo .marginal.is_ResLo
249 DemSat .marginal.is_DemSat
250 Cuts .marginal.is_Cuts /
251 dict_f 'forward simulation'
252 / is .scenario.''
253 so .opt .''
254 RES .fixed .is_RES
255 inflow .param .is_inflow
256 ACOST .level .is_ACOST
257 ALPHA .level .is_ALPHA
258 RES .level .is_RES
260 /;
261
263 Scalar start;
264 start = jnow;
265
266 option solveOpt = clear;
267 loop(j$((jnow - start)*24 < 2), // main iteration
268 is(i,s) = yes; // all realizations and trials
269 loop(revt(wp,wloop), // Backward recursion
270 tt(t) = no; tt(t)$wt(wloop,t) = yes;
271 ww(w) = no; ww(wloop) = yes;
272
273 option clear = is_RES, clear = is_inflow;
274
275 is_RES(i,s,t)$last(wloop-1,t) = i_res(i,t);
276 is_inflow(i,s,t)$wt(wloop,t) = sw_Inflow(wloop,s)/card(h);
277
278 solve waterSDDP min ACOST using lp scenario dict_b;
279
280 cont_m(j,i,wloop) = sum((s,first(wloop,t)), is_Cont(i,s,t))/card(
s);
281 delta(j,i,wloop-1) = [sum{(s,wt(wloop,t)), is_Cont(i,s,t)*is_inf
low(i,s,t)
282 + sum(ft$capacity(t,ft), is_PlantCap(i,s,t
,ft)*(-capacity(t,ft)))
283 + is_ResUp(i,s,t)*(-resmax(t)) + is_ResLo(
i,s,t)*resmin(t)
284 + is_DemSat(i,s,t)*(demand(t) - exchange(t
))}
285 + sum((s,jj,ip,t)$last(wloop,t),
286 is_cuts(i,s,jj,ip,wloop,t)*delta(jj,ip,w
loop))$(ord(wloop) < card(w))]/card(s);
287 jj(j) = yes; // we want the cuts from the stage w+1
288 ); //end of backward recursion
289
290 loop(wloop, // Forward simulation
291 tt(t) = no; tt(t)$wt(wloop,t) = yes;
292 ww(w) = no; ww(wloop) = yes;
293 if(sameas('w1',wloop),
294 RES.fx(t)$last(wloop--1,t) = resmin(t);
295 // we are using the original inflow first stage, not one of the 1
2 scenarios
296 solve waterSDDP min ACOST using lp;
297
BATINCLUDE \\ug.kth.se\dfs\home\p\v\pvshinde\appdata\xp.V2\Documents\gamsdir\pro
jdir\storeSol.gms
300 conv(j,'lo') = ACOST.l; //store the first stage cost as the low
er bound
301 zt(j,i) = ACOST.l - sum(last(ww,tt(t)), ALPHA.l(t+1));
302 i_res(i,tt)$last(wloop,tt) = RES.l(tt);
303 RES.lo(t)$last(wloop--1,t) = 0;
304 RES.up(t)$last(wloop--1,t) = inf; // free out fixings
305 else
306 is(i,s) = jwis(j,wloop,i,s); // realizations and trials sample
307 option clear = is_RES, clear = is_inflow;
308 // fix the reservoir levels of the previous stage when going forw
ard
309 is_RES(is(i,s),t)$last(wloop-1,t) = i_res(i,t);
310 is_inflow(is(i,s),t)$wt(wloop,t) = sw_Inflow(wloop,s)/card(h);
311 solve waterSDDP min ACOST using lp scenario dict_f;
312
BATINCLUDE \\ug.kth.se\dfs\home\p\v\pvshinde\appdata\xp.V2\Documents\gamsdir\pro
jdir\storeSol.gms
315 i_res(i,tt)$last(wloop,tt) = sum(is(i,s), is_RES(i,s,tt));
316 zt(j,i) = zt(j,i) + sum(is(i,s), is_ACOST(i,s) - sum(last(ww,tt(t
)), is_ALPHA(i,s,t+1)));
317 ); // end of if
318 ); // end for forward simulation
319 conv(j,'up') = sum(i, zt(j,i))/card(i);
320
* Check convergence (original Pereira and Pinto (1991) criterion)
conv(j,'sigma') = sqrt[sum(n, sqr(conv(j,'up') - zt(j,n)))/sqr(card(n))
];
break$(conv(j,'lo') > (conv(j,'up') - conv(j,'sigma')) and
conv(j,'lo') < (conv(j,'up') + conv(j,'sigma')));
327
328 * Convergence criterion of SDDP presented by Shapiro (2009) - remark1
329 conv(j,'sigma') = sqrt[sum(i, sqr(conv(j,'up') - zt(j,i)))/(card(i) - 1
)];
330 break$(abs[conv(j,'lo') - (conv(j,'up')+1.96*conv(j,'sigma')/sqrt(card(
**** $140
i)))] < 0.01*abs(conv(j,'lo')));
$36
**** LINE 340 INPUT \\ug.kth.se\dfs\home\p\v\pvshinde\appdata\xp.V2\Doc
uments\gamsdir\projdir\SDDP_example1.gms
331 conv(j,'time [min]') = (jnow - start)*24*60;
332 );
GAMS 24.7.3 r58181 Released Jul 11, 2016 WEX-WEI x86 64bit/MS Windows 04/22/20 09:03:31 Page 3
SDDP Algorithm
Error Messages


36 '=' or '..' or ':=' or '$=' operator expected
rest of statement ignored
140 Unknown symbol

**** 2 ERROR(S) 0 WARNING(S)
GAMS 24.7.3 r58181 Released Jul 11, 2016 WEX-WEI x86 64bit/MS Windows 04/22/20 09:03:31 Page 4
SDDP Algorithm
Include File Summary


COMPILATION TIME = 3.219 SECONDS 5 MB 24.7.3 r58181 WEX-WEI


**** USER ERROR(S) ENCOUNTERED

[bussieck]

break (https://www.gams.com/beta/docs/UG_FlowC ... trol_Break) has been introduced in 24.8 (https://www.gams.com/31/docs/RN_248.htm ... ue%20break). You use GAMS 24.7.3. The sddp example has been introduced in 23.7. So I suggest you take the sddp example from your distribution, not from the web site. The web site features models that can be run with the current version of GAMS. We might use features/language elements that are not available in older distributions. You might want to consider to update the software.

-Michael