$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: https://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 https://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: symbols: all $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; );