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.
No comments:
Post a Comment