Tuesday, November 11, 2014

Progressive Hedging (1)

In Pallson and Ravn (http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=6738) two small, basic examples of the progressive hedging algorithm are given. It is always a good exercise to see if we can reproduce the results. The first example is actually not very difficult to reproduce:

Algorithm:

image

Problem:

image

 

$ontext

 
Olafur P. Pallson, Hans F. Ravn
 
Scenario Analysis and the Progressive Hedging Algorithm
   
- Simple Numerical Examples
 
The Technical University of Denmark, 1993


$offtext

set
   iter
'iteration number' /it0*it9/
   s   
'scenarios' /s1*s2/
;

table ds(s,*) 'single stochastic parameter'
     
value   prob
 
s1    5     0.6
 
s2    2     0.4
;

parameter
   Xs(s)
   Ws(s) 
'price'
   report(iter,*,*)
'table 1: results from iterations'
;
Ws(s) = 0;

scalars
   iterno
   continue
/1/
   d  
/0/
   xhat
'conditional expectation' /0/
   r 
'penalty parameter'/0/
   w 
/0/
;


variables
   x 
'for submodel'
   z 
'for submodel'
;

equations
   obj
;


obj.. z =e= sqr(x-d) + w*x + 0.5*r*sqr(x-xhat);
x.lo=3;
x.up=6;
model m /obj/;

m.solprint = 2;
m.solvelink = %Solvelink.LoadLibrary%;

loop(iter$continue,

  iterno =
ord(iter);
 
display '***** iteration ******',iterno;

* r = 0 in first iteration (it0)
  r = 1$(iterno>1);

* step 1
* solve each scenario
 
loop(s,
     d = ds(s,
'value');
     w = ws(s);
    
solve m minimizing z using nlp;
     Xs(s) = x.l;
  );
 
display Xs;

* step 2
* Calculate policy Xhat
  Xhat =
sum(s, Xs(s)*ds(s,'prob'));
 
display Xhat;

* step 3
* Update W
  Ws(s) = Ws(s) + r*(xs(s) - xhat);
 
display Xhat,Ws;


  report(iter,
'x',s) = Xs(s);
  report(iter,
'xhat','-') = xhat;
  report(iter,
'w',s) = Ws(s);
  report(iter,
'z','-') = sum(s,ds(s,'prob')*sqr(Xs(s)-ds(s,'value')));
  report(iter,
'zhat','-') = sum(s,ds(s,'prob')*sqr(xhat-ds(s,'value')));

);

option report:2:1:2;
display report;



*
* comparison
* this should give the same solution
*

equation obj2;
obj2.. z =e=
sum(s,ds(s,'prob')*sqr(x-ds(s,'value')));
model m2 /obj2/;
solve m2 minimizing z using nlp;
display x.l;

 

Results in paper:

image

Results from GAMS model:

image

xhat converges to the optimal solution xhat=3.8 quite quickly.