GAMS optimizing assignments

The following fragment shows an unexpected performance problem. I thought all these assignments should be fast (i.e. only over the nonzero elements).

   1  set i /i1*i300/;
   2  alias (i,j,k);
   3   
   4  parameter p(i,j,k);
   5  p(i,i,i) = 3.14;
   6   
   7  parameter p1(i,j,k),p2(i,j,k);
   8   
   9  p1(i,j,k) = -p(i,j,k);
  10  p2(i,j,k)$p(i,j,k) = -p(i,j,k);
  
----      1 ExecInit                 0.000     0.000 SECS      3 Mb 
----      5 Assignment p             0.000     0.000 SECS      4 Mb    300
----      9 Assignment p1            2.683     2.683 SECS      4 Mb    300
----     10 Assignment p2            0.000     2.683 SECS      4 Mb    300

We see here the assignment to p1 is slow (GAMS runs through the carthesian product i x j x k) while the assignment to p2 is fast. I would expect the assignment to p1 to be fast as well.

An assignment of the form
p1(i,j,k) = (-1)*p(i,j,k)
would be optimized by GAMS while
p1(i,j,k) = -1*p(i,j,k)
would be slow again.

The profiler can help you detect these situations and optimize the assignments manually like I did for p2.

Hessian via PathNLP

Same model as Hessian via Convert but now using PathNLP.

$ontext
 
 
    Variance estimator for an MLE
 
    Erwin Kalvelagen July 2008
 
    Reference:
       William H. Greene, "Econometric Analysis", 5th edition, 2003
 
$offtext
 
set i 'cases' /i1*i20/;
 
table data(i,*)
 
                Y         E
     i1         20.5      12
     i2         31.5      16
     i3         47.7      18
     i4         26.2      16
     i5         44.0      12
     i6         8.28      12
     i7         30.8      16
     i8         17.2      12
     i9         19.9      10
    i10         9.96      12
    i11         55.8      16
    i12         25.2      20
    i13         29.0      12
    i14         85.5      16
    i15         15.1      10
    i16         28.5      18
    i17         21.4      16
    i18         17.7      20
    i19         6.42      12
    i20         84.9      16
;
 
parameter y(i),x(i);
y(i) = data(i,'y');
x(i) = data(i,'E');
 
scalar Var2Greene 'Variance estimate number two from Greene' / 46.16337 /;
 
variable b,L;
equation loglik;
 
loglik.. L =e= sum(i,log(b+x(i))) + sum(i,y(i)/(b+x(i)));
 
 
*
* estimation via pathnlp
*
 
$onecho > pathnlp.opt
hessian_gdx 1
$offecho
 
option nlp=pathnlp;
model estimation /loglik/;
estimation.optfile=1;
solve estimation minimizing L using nlp;
display "estimate",b.l;
 
 
*
* get hessian
*
scalar d2l__b_b 'hessian';
execute_load "hesData.gdx",d2l__b_b;
display d2l__b_b;
 
*
* inverting a 1x1 matrix is not very difficult
*
scalar v 'variance';
v = 1/d2l__b_b;
display v,var2greene;
 
abort$(abs(v-var2greene)>1.0e-5) "Accuracy check failed";

PathNLP prints all kind of debug output. This seems not completely polished yet.

Hessian via Convert

GAMS has new facilities to retrieve the Hessian matrix (matrix of second derivatives). Here we use Convert to produce a gdx file with the Hessian which we use to estimate the Variance of a Maximum Likelihood Estimation result.

$ontext
 
 
    Variance estimator for an MLE
 
    Erwin Kalvelagen July 2008
 
    Reference:
       William H. Greene, "Econometric Analysis", 5th edition, 2003
  
$offtext
  
set i 'cases' /i1*i20/;
  
table data(i,*)
 
                Y         E
     i1         20.5      12
     i2         31.5      16
     i3         47.7      18
     i4         26.2      16
     i5         44.0      12
     i6         8.28      12
     i7         30.8      16
     i8         17.2      12
     i9         19.9      10
    i10         9.96      12
    i11         55.8      16
    i12         25.2      20
    i13         29.0      12
    i14         85.5      16
    i15         15.1      10
    i16         28.5      18
    i17         21.4      16
    i18         17.7      20
    i19         6.42      12
    i20         84.9      16
;
 
parameter y(i),x(i);
y(i) = data(i,'y');
x(i) = data(i,'E');
 
scalar Var2Greene 'Variance estimate number two from Greene' / 46.16337 /;
 
variable b,L;
equation loglik;
 
loglik.. L =e= sum(i,log(b+x(i))) + sum(i,y(i)/(b+x(i)));
 
 
*
* estimation
*
model estimation /loglik/;
solve estimation minimizing L using nlp;
display "estimate",b.l;
 
 
*
* get hessian
*
option nlp=convert;
$onecho > convert.opt;
hessian
$offecho
estimation.optfile=1;
solve estimation minimizing L using nlp;
 
*
* gams cannot add elements at runtime so we declare the necessary elements here
*
set dummy /e1,x1/;
 
parameter h(*,*,*) '-hessian';
execute_load "hessian.gdx",h;
display h;
 
*
* inverting a 1x1 matrix is not very difficult
*
scalar v 'variance';
v = -1/h('e1','x1','x1');
display v,var2greene;
 
abort$(abs(v-var2greene)>1.0e-5) "Accuracy check failed";

Note: if you get a message

*** Msg=Error loading library c:\progra~1\gams\g2d: The specified module could not be found.

then copy gd.dll from the GAMS system directory to c:\progra~1\gams. This bug is already fixed.

GDX question

> I would like to make a gdx dump of all the sets/paramters/scalars of
> my model, but I do not want any equations or variables in there at
> all.

$version 227
 
  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 ;
 
 
* generate include file execute_unload statement
* that write all symbols except variables and
* equations
$onecho > f.awk
BEGIN {
print "execute_unload 'mygdx'" > "temp1.inc"
}
/^ *[0-9]+ / {
  if ($4=="Equ") next;
  if ($4=="Var") next;
  print "," $2 > "temp1.inc"
}
END {
print ";" > "temp1.inc"
}
$offecho
$gdxout temp1
$unload
$gdxout
$call 'gdxdump temp1 symbols > temp1.txt'
$call 'awk -f f.awk temp1.txt'
$include temp1.inc
 
* mygdx.gdx will not contain any variables or equations
execute 'gdxdump mygdx symbols';

Basically at compile time we dump a gdx file temp1.gdx that has all symbols. Then we use an AWK scripts (still at GAMS compile time) to generate an include file temp1.inc that contains an execute_unload statement with all non var/equ symbols to mygdx.gdx. Then we $include this file. Finally for checking we call gdxdump (at execution time) to show all symbols in mygdx.gdx.

gdx2xls does not work with 22.7

> Erwin. When tried to create xls file from gdx, I got an error message (see log file).
>*** Error at line 7379: Execute has nonzero return code RC=2

The GDX API has been changed between 22.5 and 22.6. You hare using 22.7 and your version of gdx2xls for <= 22.5. A work around is:

scalar s /3.14/;
execute_unload "s.gdx",s;
execute '=gdxcopy -v6c s.gdx oldgdx';
execute '=gdx2xls oldgdx\s.gdx';
execute '=shellexecute oldgdx\s.xls';.

This should show:

--- Job Untitled_39.gms Start 07/24/08 18:17:48
GAMS Rev 227  Copyright (C) 1987-2008 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen                               G070509/0001CE-WIN
          GAMS Development Corporation                               DC4572
--- Starting compilation
--- Untitled_39.gms(9) 2 Mb
--- Starting execution: elapsed 0:00:00.003
--- Untitled_39.gms(7) 3 Mb
Writing files to: .\.\oldgdx
Converting file: .\s.gdx
--- Untitled_39.gms(8) 3 Mb
GDX2XLS Version 3.0
Reading GDX file oldgdx\s.gdx
Reading 1 symbols, sorting: 0.00 seconds
Reading 0 UELs: 0.00 seconds
Opening C:\projects\test\oldgdx\s.XLS with Excel: 0.31 seconds
           s. To scalar sheet
Total elapsed time: 0.72 seconds
--- Untitled_39.gms(9) 3 Mb
Shell Execute    May  1, 2008 22.7.2 WIN 3727.4799 VIS x86/MS Windows         
--- Untitled_39.gms(9) 3 Mb
*** Status: Normal completion
--- Job Untitled_39.gms Stop 07/24/08 18:19:49 elapsed 0:02:00.857

Update: unfortunately this worked for me, but not for the user. I don't understand why. Contact support@gams.com and ask for a 22.7 version of gdx2xls.exe.

gdxcopy puzzle

This does not look right:

c:\projects\test>gdxcopy -v6c save\s.gdx oldgdx
Writing files to: save\save\oldgdx
Converting file: save\s.gdx

Instead of writing to .\oldgdx or save\oldgdx it goes to save\save\oldgdx.

Update: this is fixed in 22.8 beta. Minor issue: the log messages have disappeared in the 22.8 version.

c:\projects\test>gdxcopy -v6c save\s.gdx oldgdx
 
c:\projects\test>

I would prefer to see what gdxcopy is doing like in the previous version.

Modeling question

>Among the constraints in the formulation my team have worked on,
>some of them have special properties like:
>x1>=a1
>x1+x2>=a1+a2
>x1+x2+x3>=a1+a2+a3
>...
>...
>x1+x2+x3+...+x39>=a1+a2+a3+...+a39
>
>Here a[i] are nonnegative parameters, so RHSs are non-decreasing.
>Besides these constraints, we deal with tons of constraints, so we are not looking
>for any algorithm or method, but looking for smart conversion of these constraints.
>My question is "Are they any ways to formulate this more intellingently or effectively
>to reduce computational effort?


parameter b(i);
alias (i,j);
b(i) = sum(j$(ord(j)<=ord(i)), a(j));

variable y(i);

ydef(i).. y(i) =e= y(i-1) + x(i);
y.lo(i) = b(i);

This requires more variables but many fewer nonzeroes, so likely more efficient (sparser).

Excel GUI for GIS and GAMS

This is work in progress: Excel as front-end for a GAMS Model, including a GIS component. Click to enlarge.

Modeling question

>I need to change the following non linear constraint to a linear one .
>the constraints are:
>
>for all k = 0 ,1, 2,3 ,4,...... 10
>
>k+1 = sum over j ( a[k][j] * d[j] )
>
>j = 1,2,3,4
>
>for all k, j a[k][j] is an integer <= 2 and d[j] is an integer <= 60. >
>both a[k][j] and d[j] are variables , making it a non linear
>constraint. Is there a method i can change it to a linear expression.
>Cplex gives a error with this constraint.

These question can only be answered more definitely by looking at the complete model, but let me give it a shot anyway.

The a(k,j) variables can only assume the values 0,1,2. So we can write this as a sum of 2 binary variables (and relax a to be continuous):

0 <= a(k,j) <= 2 continuous
b1(k,j), b2(k,j) binary
a(k,j) = b1(k,j)+2*b2(k,j);

Now we can express z(k,j) = a(k,j)*d(j) as

M1=20
z1(k,j) ≤ M1 · b1(k,j)
z1(k,j) ≤ d(j)
z1(k,j) ≥ d(j) − M1 · (1 − b1(k,j))
M2=40
z2(k,j) ≤ M2 · b2(k,j)
z2(k,j) ≤ 2*d(j)
z2(k,j) ≥ 2*d(j) − M2 · (1 − b2(k,j))
k + 1 = sum(j, z1(k,j)+z2(k,j));

See also: http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-continuous-and-binary.html.

Invert.exe raw mode gdx

Using raw mode to write the inverse (instead of string mode) yields:

02.416 seconds (n=1000)

The GAMS version should get similar performance (may be even faster as I did not use any compiler flags to do better optimization).