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:

Problem:

 \$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:

Results from GAMS model:

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