Progressive Hedging (2)

The second example in Pallson and Ravn (http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=6738) is a little bit more challenging.

We have:

with the following tree:

A GAMS representation of the progressive hedging algorithm could be:

$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*it20/    n       'nodes'     /I,II,III,IV, s1*s6/    s(n)    'scenarios: leaf nodes' /s1*s6/    t       'levels: dimension of x' /t1*t3/    dp(n)   'decision points in tree' /I,II,III,IV/    tt(t)   "subset of t's used in current scenario" ; alias(n,nn,nnn); table tree(n,nn,*)         value  prob I.II      5    0.3 II.s1     6    0.2 II.s2     5    0.5 II.s3     2    0.3 I.III     2    0.7 III.s4    4    0.6 III.IV    7    0.4 IV.s5     6    0.1 IV.s6     4    0.9 ; display tree; * * do some tree manipulations here * set paths(s,n,nn) 'links on path to leaf node'; paths(s,n,s)$tree(n,s,'prob')=yes; scalar cnt; * repeat augment paths until no more change repeat     cnt = card(paths);     paths(s,nnn,n)$(sum(nn,paths(s,n,nn)) and tree(nnn,n,'prob')) = yes; until cnt=card(paths); option paths:0:0:1; display paths; set top(n) 'top of tree'; top(n)$(sum((s,nn)$paths(s,nn,n),1)=0) = yes; display top; parameter level(n) 'level of each node in tree'; level(top)=1; repeat    cnt = card(level);    loop(paths(s,n,nn)$level(n),       level(nn) = level(n)+1;    ); until cnt=card(level); display level; alias(t,tdown,tup); parameters    sprob0(s,t) 'single probabilities'    sprob1(s,t) 'accumulated probabilities'    sval(s,t) 'values' ; sval(s,t) = sum(paths(s,n,nn)$(level(n)=ord(t)), tree(n,nn,'value')); display sval; sprob0(s,t) = sum(paths(s,n,nn)$(level(n)=ord(t)), tree(n,nn,'prob')); display sprob0; * accumulate probabilities backwards to get conditional probabilities sprob1(s,t)$sprob0(s,t) = prod(tdown$(ord(tdown)>=ord(t) and sprob0(s,tdown)), sprob0(s,tdown)); display "table 4a",sprob1; set stdp(s,t,dp) 'mapping to decision points'; stdp(s,t,dp)$sum(paths(s,dp,nn)$(level(dp)=ord(t)), 1) = yes; option stdp:0:1:2; display stdp; set tdp(t,dp) 'extracted from stdp for a given scenario'; parameter    Xs(s,t)    Ws(s,t)  'price'    W(t)    d(t)    xhat(dp) 'conditional expectation'    report(iter,*,*) 'table 1: results from iterations'    x0(t) /t1 4/ ; Ws(s,t) = 0; w(t) = 0; xhat(dp) = 0; scalars    iterno    continue /1/    r  'penalty parameter'/0/ ; variables    x(t)  'for submodel'    z  'for submodel' ; equations    obj ; obj.. z =e= sum(tt(t),0.5*sqr(x(t)-x(t-1)-x0(t))+sqr(x(t)-d(t)))      + sum(tt(t),w(t)*x(t)) + 0.5*r*sum(tt(t),sqr(x(t)-sum(tdp(t,dp),xhat(dp)))); x.lo(t)=3; x.up(t)=6; model m /obj/; * reduce output and increase speed in loop m.solprint = 2; m.solvelink = %Solvelink.LoadLibrary%; option xhat:2; display xhat; set report_iter(iter) /it0*it5,it10,it15,it20/; parameter report_xhat(iter,dp) 'table 7'; loop(iter$continue,   iterno = ord(iter);   display '***** iteration ******',iterno; * r = 0 in first iteration (it0)   r = 1$(iterno>1); * step 1 * solve each scenario   Xs(s,t)=0;   loop(s,      d(t) = sval(s,t);      tdp(t,dp) = stdp(s,t,dp);      tt(t) = sum(tdp(t,dp),1);      w(t) = ws(s,t);      solve m minimizing z using nlp;      Xs(s,tt(t)) = x.l(t);   );   display Xs; * step 2 * Calculate policy Xhat   xhat(dp) = sum(stdp(s,t,dp), Xs(s,t)*sprob1(s,t));   display Xhat; * step 3 * Update W   Ws(s,t) = Ws(s,t) + r*(xs(s,t) - sum(stdp(s,t,dp),xhat(dp)));   display Ws;   report_xhat(iter,dp)$report_iter(iter) = xhat(dp); ); display report_xhat;

Much of the “ugly” code is dealing with the tree. We need to manipulate it to get the following data:

This is reproduced in GAMS:

The PHA algorithm itself is largely the same as for the smaller first problem shown here: http://yetanothermathprogrammingconsultant.blogspot.com/2014/11/progressive-hedging-1.html Interestingly we have some differences in the behavior of the algorithm. The paper shows the following iteration results:

while I get:

Apart from iteration 0, we are not doing the same thing. But looking at the optimal solution for xhat:

I am actually a bit closer to this optimal solution.

We cannot see what the differences are between our calculations and the ones performed in the paper. It would have been easier if the code the authors used to get to their results would be available as an appendix.

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.

Process GAMS code

For a documentation tool, I investigated how we can use compiled GAMS code to understand data manipulation statements and model equations. I have seen attempts in the past where GAMS source was parsed directly. This is is a very difficult task. The GAMS language has many esoteric extensions and maintaining full compatibility is almost impossible. Indeed in the examples I looked at even small standard models from the model library were rejected.

So instead of parsing GAMS models directly, I work with the output of the GAMS compiler. This is a byte-code type of output, that is not very suited for further processing (other than execution). So we take this linear code and make a tree out of it. Walking the tree and generating all kind of interesting output is then relatively straightforward. Here we produce LaTeX output. E.g.

c(i) = sum(j, f(j)*d(i,j))/100;

causes the GAMS compiler to generate:

INSTRUCTION DUMP FROM 1 TO 401

  LOC    NUM INSTRUCT   SUB        FLD IDENT

  114      0 UnitBeg       0         17
  115     14 DefBeg        0        137 c
  116      7 CntrBeg       0        133 i
  117     12 Index         0          1
  118     55 EndLhs        0          0
  119     29 SumBeg        0          0
  120      7 CntrBeg       0        134 j
  121     12 Index         0          2
  122      2 Push          0        139 f
  123     12 Index         0          1
  124     12 Index         0          2
  125      2 Push          0        138 d
  126      4 MultOp        0          0
  127     30 SumOp         0          0
  128      8 CntrEnd       0          0
  129     21 Immed         0          5   1.00000000000000E+002
  130      4 MultOp        1          0
  131     15 DefOp         0        137 c
  132      8 CntrEnd       0          1
  133     16 DefEnd        0        137 c
  134      1 UnitEnd       0          0

A tree version of this can look like:

Finally some output:

A more elaborate equation:

armington(i)..
  x(i) =e= ac(i)*(delta(i)*m(i)**(-rho(i)) + (1-delta(i))*xxd(i)**(-rho(i)))**(-1/rho(i)) ;

This gives a large tree:

If we also use colors to indicate if a symbol is a variable (red) or a parameter (blue) we can generate:

GAMS –> R

Handling of inventory

In this book and the accompanying web site we see some unusual ways to model inventory. In the first chapters, inventory is not dealt with explicitly and therefore we look at equations such as:

I.e. we add up all production up to now, add initial inventory and subtract all demand up to now (either final demand by customers or intermediate demand because of production of other products).

This is in general not how we would like to write such equations. It is much better to explicitly introduce inventory variables I_i,t≥0 and write:

I.e. inventory = previous inventory + production – demand. We add variables but the model becomes much sparser (fewer nonzero elements).

Even when they add inventory explicitly in models later on, they still do it wrong. They write:

This is not how we want it, as we repeat many expressions over different equations. A better formulation is:

There are different ways to write down the initial inventory.  Here are some possibilities:

I have a preference for the last formulation.

MIQP vs MIP

On my experience many MIQP models can and should be reformulated into straight linear MIP models. The reason is threefold: MIQP solvers are not as reliable as MIP solvers (they can experience numerical difficulties much more often), they can not solve many QPs that are not convex and the performance of a corresponding MIP model can be much better.

An example of this was discussed in https://groups.google.com/forum/#!topic/gurobi/t-N23d24uuw.

With some larger data sets (randomly generated), I saw the following:

Scip (MIQP formulation) was stopped when exceeding 1000 seconds.

This is a non-trivial reformulation, so the presolvers are not (yet) able to replicate this.

PS. In general I don’t think it is a good idea to skimp on variables and equations as is suggested in the discussion. Often we see poorly formulated models as a result of trying to save on variables and constraints, resulting in models that are much denser and actually more difficult to solve. From the above numbers we can see my liberal use of variables and equations in the MIP formulation does not really harm and actually this formulation is light years ahead in terms of performance and reliability compared to the original MIQP formulation.

GAMS to Python

When we need to export data from GAMS to Python for further processing, I often use SQLite as intermediate storage. It is very easy to load data from SQLite into a Pandas dataframe:

 

The text_factory thing may help with Unicode issues.

Compare this to reading the same data into R:

For more info about doing data stuff in Python see:

Data manipulation in GAMS, R or SQL

When moving data from GAMS to R via SQLite we have three possible spots where we can apply some data manipulation:

Suppose we want to plot some data but need to do:

Here is the code in GAMS:

I could do the same thing in R, but the code is slightly more messy:

The final attempt is in SQL, which is not too bad either (but with nested queries more complicated than the GAMS version):

 

In all cases we did not use any loops.

Generating Box Plots from GAMS

Of course using R’s ggplot2 package. Used gdx2sqlite to pass on data from GAMS to R. See also: R Maps from GAMS,

Newer version:

Experiment: Database Diagram of GAMS model

It is always helpful to look at your model in a different way. Not sure if this is a very useful tool in practice. But it may help when passing on data from a GAMS model to a DBA as a documentation tool.

Sets        i   canning plants   / seattle, san-diego /        j   markets          / new-york, chicago, topeka / ;   Parameters        a(i)  capacity of plant i in cases          /    seattle     350               san-diego   600  /        b(j)  demand at market j in cases          /    new-york    325               chicago     300               topeka      275  / ;   Table d(i,j)  distance in thousands of miles                     new-york       chicago      topeka       seattle          2.5           1.7          1.8       san-diego        2.5           1.8          1.4  ;   Scalar f  freight in dollars per case per thousand miles  /90/ ;   Parameter c(i,j)  transport cost in thousands of dollars per case ;             c(i,j) = f * d(i,j) / 1000 ;   Variables        x(i,j)  shipment quantities in cases        z       total transportation costs in thousands of dollars ;   Positive Variable x ;   Equations        cost        define objective function        supply(i)   observe supply limit at plant i        demand(j)   satisfy demand at market j ;   cost ..        z  =e=  sum((i,j), c(i,j)*x(i,j)) ;   supply(i) ..   sum(j, x(i,j))  =l=  a(i) ;   demand(j) ..   sum(i, x(i,j))  =g=  b(j) ;   Model transport /all/ ;   Solve transport using lp minimizing z ;   Display x.l, x.m ;