Unit Loss Function

A standard normal unit loss function can be described as:

  

where f(x) is the standard normal density function and F(x) is the standard normal CDF. See: http://planetmath.org/unitnormallossfunction.

Often we want to solve UNL(x)=c for x.

GAMS Version

$ontext    Standard Normal Unit Loss Function (UNL)    solve UNL(x)=c    for a given c    Method: solve f(x) - x * (1-F(x)) = c    where f is standard normal density function and F is the standard    normal CDF. $offtext scalar c /0.0331/; variable x; equation e; e..  1/sqrt(2*pi)*exp(-0.5*sqr(x)) - x * (1-errorf(x)) =e= c; * initial value x.l = 1; model m /e/; solve m using cns; display x.l;

This solves with CONOPT very quickly, we get:

Note: the GAMS math expressions can be inspected with a new tool. This is a nice test example:

Wolfram Alpha

The input is easy:

The output is unfortunately:

Mathematica

The above approach is not really suited for Mathematica. This looks better:

R

Super easy input, due to predefined functions for normal pdf and cdf. We give an interval of [0,100] to search:

Excel

Its root finder is called Goal Seek.

It is a little bit sloppy by default. We can make Goal Seek more precise by multiplying the UNC-c cell by a 1000:

  

Note: the tolerance can also be set in the Excel options (under formulas, Maximum Change). Looks like Excel is using an absolute tolerance on Δx and Δf (i.e. it could stop prematurely if the algorithm chooses a very small Δx).

Python

Quite nice:

Quartiles in R

In this Wikipedia article it looks like the R functions summary and fivepoint use the same algorithm to calculate q1 and q3:

However that is not the case, as can be seen here (notice the difference in q1 and q3):

Of course quartiles are not uniquely defined. In fact the function quantile offers nine(!) different ways to calculate them: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/quantile.html. The function summary actually calls quantile using its default algorithm type 7.

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.