sddp.gms : Multi-stage Stochastic Water Reservoir Model solved with SDDP

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:
    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;
);