Description
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.
Large Model of Type : LP
Category : GAMS Model library
Main file : sddp.gms includes : water.gdx water.xlsx
$title Multi-Stage Stochastic Water Reservoir Model solved with SDDP (SDDP,SEQ=357)
$onText
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, GUSS,
water resource management
$offText
$eolCom //
* Only test fast LP
$ifI %gams.lp%==cplex $goTo continue
$ifI %gams.lp%==xpress $goTo continue
$ifI %gams.lp%==gurobi $goTo continue
$log SDDP only tested for selected LP solvers
$exit
$label continue
Set
w 'weeks' / w1*w52 /
t 'hours of the year' / t1*t8736 /
ft 'types of fuel available (plants)' / Hydro, HardCoal, Nuclear /
s 'scenarios' / s1*s12 /;
Parameter
demand(t) 'power demand MW'
exchange(t) 'cross border flow MW'
wCapacity(w,ft) 'capacity of plant MW'
wPrices(w,*) 'coal and CO2 price EUR per ton'
wInflow(w) 'inflows into the reservoir MWh per week'
sw_Inflow(w,s) 'inflow realizations MWh per week'
resmax(t) 'max reservoir level MW' / t1*t8735 106.2e6, t8736 87e6 /
resmin(t) 'min reservoir level MW' / t1*t8735 10e6, t8736 87e6 /;
$onEmbeddedCode Connect:
- ExcelReader:
file: water.xlsx
symbols:
- {name: demand, rowDimension: 1, columnDimension: 0, range: t_Demand!A1}
- {name: exchange, rowDimension: 1, columnDimension: 0, range: t_Flow!A1}
- {name: wCapacity, rowDimension: 1, columnDimension: 1, range: w_Capacity!A2}
- {name: wPrices, rowDimension: 1, columnDimension: 1, range: w_Prices!A2}
- {name: wInflow, rowDimension: 1, columnDimension: 0, range: w_Inflow!A1}
- {name: sw_Inflow, rowDimension: 1, columnDimension: 1, range: sw_Inflow!A2}
- GAMSWriter:
writeAll: True
$offEmbeddedCode
* Week hour mapping
Set
h 'hours of one week' / h1*h168 /
wt(w,t) 'map weeks and hours'
last(w,t) 'last hour t of each week w'
first(w,t) 'first hour t of each week w';
loop(w,
loop((h,t)$sameas('t1',t), wt(w,t + ((ord(w) - 1)*card(h) + ord(h) - 1)) = yes;);
loop(t$sameas('t1',t),
first(w,t+((ord(w) - 1)*card(h))) = yes;
last(w,t+(ord(w)*card(h) - 1)) = yes;
);
);
Parameter
capacity(t,ft) 'capacity of plant MW'
gencost(t,ft) 'generating cost EUR per MW' / #t.Hydro 0, #t.Nuclear 15 /
inflow(t) 'inflows into the reservoir MW';
inflow(t) = sum(wt(w,t), wInflow(w)/card(h));
capacity(t,ft) = sum(wt(w,t), wCapacity(w,ft));
gencost(t,'HardCoal') = sum(wt(w,t), (wPrices(w,'HardCoal') + 2.361*wPrices(w,'CO2')))/(6.98*0.39);
Set
tt(t) 'control set for hours'
ww(w) 'control set for weeks';
Positive Variable
GAP(t) 'gap generation in hour t MW'
RES(t) 'reservoir energy (level) at end of t MW'
SPILL(t) 'spillage in hour t MW'
X(t,ft) 'generated power in hour t by plant with fuel type MW'
SLACKUP(t) 'slack variable for the upper reservoir level equation MW'
SLACKLO(t) 'slack variable for the lower reservoir level equation MW';
Variable COST 'total generation cost EUR';
Equation
Cont(t) 'hydraulic continuity equation'
DemSat(t) 'demand satisfaction'
ResUp(t) 'maximum reservoir level'
ResLo(t) 'minimum reservoir level'
PlantCap(t,ft) 'plant power generation capacity'
Obj 'objective function';
Cont(tt(t))..
RES(t) - RES(t--1) + X(t,'Hydro') + SPILL(t) =e= inflow(t);
DemSat(tt(t))..
sum(ft$capacity(t,ft), X(t,ft)) + GAP(t) =g= demand(t) - exchange(t);
ResUp(tt(t))..
-RES(t) + SLACKUP(t) =g= -resmax(t);
ResLo(tt(t))..
RES(t) + SLACKLO(t) =g= resmin(t);
PlantCap(tt(t),ft)$capacity(t,ft)..
-X(t,ft) =g= -capacity(t,ft);
Obj.. COST =e= sum((tt(t),ft)$capacity(t,ft), gencost(t,ft)*X(t,ft))
+ sum(tt(t), 1000*GAP(t) + 10e6*(SLACKUP(t) + SLACKLO(t)));
Model water / all /;
* Deterministic model
tt(t) = yes;
$if set solvedet solve water using lp min Cost;
$title SDDP Algorithm
* Scenario data
$if not set jmax $set jmax 20
$if not set imax $set imax 5
Set
i 'trial solutions' / i1*i%imax% /
j 'iteration index' / j1*j%jmax% /
jj(j) 'dynamic j';
Parameter
cont_m(j,i,w) 'dual variables associated with the mass balance constraint'
delta(j,i,w) 'RHS of the Benders cuts';
cont_m(j,i,w) = 0;
delta(j,i,w) = 0;
jj(j) = no;
Free Variable ACOST 'approximation of cost';
Positive Variable ALPHA(t) 'approximation of future cost function (FCF)';
Equation
Obj_Approx 'objective function for the one-stage dispatch model'
Cuts(j,i,w,t) 'renders cuts';
Obj_Approx..
ACOST =e= sum((tt(t),ft)$capacity(t,ft), gencost(t,ft)*X(t,ft))
+ sum(tt(t), 1000*GAP(t) + 10e6*(SLACKUP(t) + SLACKLO(t)))
+ sum(last(ww,tt(t)), ALPHA(t+1));
Cuts(jj,i,last(ww(w),t))$(ord(w) < card(w))..
ALPHA(t+1) - cont_m(jj,i,w+1)*RES(t) =g= delta(jj,i,w);
Model waterSDDP / water, obj_approx, cuts /;
* No superflous output and in core solving
option limRow = 0, limCol = 0, solPrint = silent, solveLink = %solveLink.loadLibrary%;
Parameter i_res(i,t) 'trial decisions';
Set jwis(j,w,i,s) 'scenario trial sample';
* Select some for the first backward recursion;
loop(i, i_res(i,t) = resmin(t) + (ord(i) - 1)/(card(i) - 1)*(resmax(t) - resmin(t)));
loop(s$sameas('s1',s), jwis(j,w,i,s + UniFormInt(0,card(s) - 1))$(not sameas('w1',w)) = yes);
* We need to iterate w from back to front
Set revt(w,w) 'backward loop for backward recursion excluding week 1';
loop(w$(ord(w) < card(w)), revt(w,w + (card(w) - 2*ord(w) + 1)) = yes);
$if not set debug $set debug 0
Parameter
conv(j,*) 'convergence parameters'
zt(j,i) 'objective for trial solution';
$if %debug% == 1 sol(*,j,i,t,*) Solution matrix
$onEchoV > storeSol.gms
$ifThen %debug% == 1
$onDotL
$if not x%1==x loop(is(i,s),
sol('res',j,i,tt,'') = %1RES(%2tt);
sol('x',j,i,tt,ft) = %1X(%2tt,ft);
sol('gap',j,i,tt,'') = %1GAP(%2tt);
sol('spill',j,i,tt,'') = %1SPILL(%2tt);
sol('slackLo',j,i,tt,'') = %1SLACKLO(%2tt);
sol('slackUp',j,i,tt,'') = %1SLACKUP(%2tt);
$if not x%1==x );
$offDotL
$endIf
$offEcho
Alias (w,wp,wloop), (i,ip);
Set is(i,s) / #i.#s /;
Parameter
is_RES(i,s,t)
is_inflow(i,s,t)
is_Cont(i,s,t)
is_PlantCap(i,s,t,ft)
is_ResUp(i,s,t)
is_ResLo(i,s,t)
is_DemSat(i,s,t)
is_Cuts(i,s,j,i,w,t)
is_ACOST(i,s)
is_ALPHA(i,s,t)
$ifThen %debug% == 1
is_X(i,s,t,ft)
is_GAP(i,s,t)
is_SPILL(i,s,t)
is_SLACKLO(i,s,t)
is_SLACKUP(i,s,t)
$endIf
so / SkipBaseCase 1, LogOption 1, UpdateType 2 /;
Set
dict_b 'backward recursion'
/ is .scenario.''
so .opt .''
RES .fixed .is_RES
inflow .param .is_inflow
Cont .marginal.is_Cont
PlantCap.marginal.is_PlantCap
ResUp .marginal.is_ResUp
ResLo .marginal.is_ResLo
DemSat .marginal.is_DemSat
Cuts .marginal.is_Cuts /
dict_f 'forward simulation'
/ is .scenario.''
so .opt .''
RES .fixed .is_RES
inflow .param .is_inflow
ACOST .level .is_ACOST
ALPHA .level .is_ALPHA
RES .level .is_RES
$ifThen %debug% == 1
X .level .is_X
GAP .level .is_GAP
SPILL .level .is_SPILL
SLACKLO .level .is_SLACKLO
SLACKUP .level .is_SLACKUP
$endIf
/;
$if not set timelim $set timelim 2 // Run for max 2 hours
Scalar start;
start = jnow;
option solveOpt = clear;
loop(j$((jnow - start)*24 < %timelim%), // main iteration
is(i,s) = yes; // all realizations and trials
loop(revt(wp,wloop), // Backward recursion
tt(t) = no; tt(t)$wt(wloop,t) = yes;
ww(w) = no; ww(wloop) = yes;
option clear = is_RES, clear = is_inflow;
is_RES(i,s,t)$last(wloop-1,t) = i_res(i,t);
is_inflow(i,s,t)$wt(wloop,t) = sw_Inflow(wloop,s)/card(h);
solve waterSDDP min ACOST using lp scenario dict_b;
cont_m(j,i,wloop) = sum((s,first(wloop,t)), is_Cont(i,s,t))/card(s);
delta(j,i,wloop-1) = [sum{(s,wt(wloop,t)), is_Cont(i,s,t)*is_inflow(i,s,t)
+ sum(ft$capacity(t,ft), is_PlantCap(i,s,t,ft)*(-capacity(t,ft)))
+ is_ResUp(i,s,t)*(-resmax(t)) + is_ResLo(i,s,t)*resmin(t)
+ is_DemSat(i,s,t)*(demand(t) - exchange(t))}
+ sum((s,jj,ip,t)$last(wloop,t),
is_cuts(i,s,jj,ip,wloop,t)*delta(jj,ip,wloop))$(ord(wloop) < card(w))]/card(s);
jj(j) = yes; // we want the cuts from the stage w+1
); //end of backward recursion
loop(wloop, // Forward simulation
tt(t) = no; tt(t)$wt(wloop,t) = yes;
ww(w) = no; ww(wloop) = yes;
if(sameas('w1',wloop),
RES.fx(t)$last(wloop--1,t) = resmin(t);
// we are using the original inflow first stage, not one of the 12 scenarios
solve waterSDDP min ACOST using lp;
$ batInclude storeSol
conv(j,'lo') = ACOST.l; //store the first stage cost as the lower bound
zt(j,i) = ACOST.l - sum(last(ww,tt(t)), ALPHA.l(t+1));
i_res(i,tt)$last(wloop,tt) = RES.l(tt);
RES.lo(t)$last(wloop--1,t) = 0;
RES.up(t)$last(wloop--1,t) = inf; // free out fixings
else
is(i,s) = jwis(j,wloop,i,s); // realizations and trials sample
option clear = is_RES, clear = is_inflow;
// fix the reservoir levels of the previous stage when going forward
is_RES(is(i,s),t)$last(wloop-1,t) = i_res(i,t);
is_inflow(is(i,s),t)$wt(wloop,t) = sw_Inflow(wloop,s)/card(h);
solve waterSDDP min ACOST using lp scenario dict_f;
$ batInclude storeSol is_ i,s,
i_res(i,tt)$last(wloop,tt) = sum(is(i,s), is_RES(i,s,tt));
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)));
); // end of if
); // end for forward simulation
conv(j,'up') = sum(i, zt(j,i))/card(i);
$onText
* 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')));
$offText
* Convergence criterion of SDDP presented by Shapiro (2009) - remark1
conv(j,'sigma') = sqrt[sum(i, sqr(conv(j,'up') - zt(j,i)))/(card(i) - 1)];
break$(abs[conv(j,'lo') - (conv(j,'up')+1.96*conv(j,'sigma')/sqrt(card(i)))] < 0.01*abs(conv(j,'lo')));
conv(j,'time [min]') = (jnow - start)*24*60;
);