ViewHar and GDX2HAR problems

Users may encounter the following problem when using ViewHar or GDX2HAR:

Cannot open input file C:\delgams\tmp2.gdx
Perhaps it is not really a GDX file !
**** Fatal GDX I/O Error = -100033

This happens when a new GAMS system (22.6 or newer) is used with an old GDX2HAR (or ViewHar) tool. The tools included in GEMPACK work with GAMS 22.5 or older. To increase the confusion, GAMS also distributes GDX2HAR and then it depends on the PATH order which version of GDX2HAR will be executed. Make sure you actually use the GAMS version of GDX2HAR and things should be ok. ViewHar is not distributed by GAMS so you need to convert the GDX to an older format first. You can use GDXCOPY for this.

The underlying problem lies in unfortunate changes in the GDX API by GAMS. It should have been possible to make these changes transparent by adding new API calls instead of changing them. In general it is poor design to change existing API calls as can be seen from this example.

How to call GAMS from Excel using VSTO?

VSTO = Visual Studio Tools for Office (see http://msdn.microsoft.com/en-us/office/aa905533.aspx). This allows you to develop .NET code for use with Office applications. Basically newer versions of Office both host VBA and .NET. VBA code is stored embedded in the .xls file, while .Net is maintained outside the spreadsheet. MS Visual Studio has excellent tools to develop, debug and deploy VSTO based applications. The VBA code presented in http://www.amsterdamoptimization.com/packaging.html can be ported to .NET/VSTO without much problem.

Nonrectangular Domains

GAMS only allows rectangular domains, ie we can use

variable flow(i,j,k);

but not

set links(i,j,k) /.../;
positive variable flow(links);

It is the responsibility of the modeler to check every occurrence of FLOW() in the model equations to make sure no FLOW(i,j,k) outside the set LINKS is referenced.

AMPL allows a better, safer approach: you can specify a set as domain:

set LINKS := ....
var flow{LINKS} >= 0;

This will check that you don't use flow outside its LINKS domain.

GAMS/SCIP license

As I understand it GAMS/COINSCIP can only be licensed to academic users. I am not a lawyer but I suspect this may violate the COIN CPL license (http://www.ibm.com/developerworks/library/os-cpl.html).

"Subject to the terms of this Agreement, each Contributor hereby grants Recipient a non-exclusive, worldwide, royalty-free copyright license to reproduce, prepare derivative works of, publicly display, publicly perform, distribute and sublicense the Contribution of such Contributor, if any, and such derivative works, in source code and object code form."

It is possible to use a different license as long as: it complies with the terms and conditions of this Agreement. As GAMS/COINSCIP is part of any GAMS distribution I would guess I qualify as a Recipient.

GAMS/IDE GDX reload

The GAMS/IDE will ask and try to update an open GDX file if the GDX file has changed. This is especially the case when the IDE gets the focus again. In some cases it will try then to reload the GDX file while GAMS is writing to this file. This results in a run-time error. The IDE does not crash: it will show a popup window with the run-time error but will continue to work otherwise.

Preferably the IDE should not reload a GDX while it is being written to.

Modeling question

>Hi all, does anybody knows how to attack a problem like the following ?
>
>Z: min (sum t(sum k(Ct,k Xt,k)))
>s.t.
>1) (sum t(sum k(Ct,k Xt,k)) GE costant
>2) (sum k(Qi,k Xt,k)) + Si,t LE Bi,t + Si,t-1 (i:1..I;t:1..T)
>3) (sum t(Xt,k)) LE 1 (t:1..T)
>4) Xt,k = (0;1)
>
>All ideas will be appreciated !!!!

A problem like this is very easily coded in a modeling system like AMPL or GAMS. If I understand the notation correctly this would look like:

binary variable x(t,k);
equations
obj
e1(i,t)
e2(k)
;
obj.. z =e=sum((t,k),C(t,k)*X(t,k));
z.lo = constant;
e1(i,t).. sum(k, Q(i,k)*X(t,k))+S(i,t) =L= B(i,t)+S(i,t-1);
e2(k).. sum(t, X(t,k)) =L= 1;

Testing for convexity

>Hi all,
>
>If I have complicated objective function which can only be evaluated
>numerically, and through theoretical analysis it's really hard to judge
>whether it is convex or not...
>
>Is there a way to test if the optimization problem is convex or not, using
>non-theoretical analysis methods?

Just by numerically evaluation or sampling it is very difficult to prove anything. There is a tool to estimate whether a problem is convex: MProbe.

If one would have a problem formulated in GAMS, the behavior of BARON can tell you if a model is convex (no branch-and-bound).

It would be a useful to have a "dummy" solver that would tell us if a problem is convex. Essentially just the first phase of BARON without actually solving the problem.

LPsolve to GAMS conversion

Some clients want to see how GAMS LP/MIP solvers perform vis-a-vis LPsolve. To this effect I created a small tool to convert LPsolve models to GAMS. Then you can run the GAMS model against any of the LP/MIP solvers available with GAMS.

An example of use is reproduced below:

D:\lpsolve\LP_SOL~1.11_>type lp.lp
max: -x1 + 2 x2;
C1: 2x1 + x2 < 5;
-4 x1 + 4 x2 <5;
   
int x2,x1;
  
  
D:\lpsolve\LP_SOL~1.11_>lp2gams lp.lp
lp2gams.exe version 1.0
Erwin Kalvelagen, erwin@amsterdamoptimization.com
   
lp_solve version 5.5.0.11
Reading lp.lp
Rows: 2   Columns: 2   Nonzero elements: 4
Writing lp.gms
Done
  
D:\lpsolve\LP_SOL~1.11_>type lp.gms
* Produced by lp2gams.exe version 1.0
* Erwin Kalvelagen, erwin@amsterdamoptimization.com
*
* Lp file: lp.lp
* Rows: 2   Columns: 2   Nonzero elements: 4
  
free variable gamsobjvar;
integer variable x1;
integer variable x2;
 
equation R0;
equation C1;
equation R2;
  
R0.. gamsobjvar =e=  - x1 + 2 * x2;
C1..  + 2 * x1 + x2 =l= 5;
R2..  - 4 * x1 + 4 * x2 =l= 5;
   
option optcr=0;
model gamsmodel /all/;
solve gamsmodel using mip maximizing gamsobjvar;
   
D:\lpsolve\LP_SOL~1.11_>
  

Median DNLP

The median (http://en.wikipedia.org/wiki/Median) can be found by minimizing sum(i, abs(x-a(i))). This is not a well-defined optimization model, and many solvers have problems with it.

An AMPL representation can be found here: http://orfe.princeton.edu/~rvdb/ampl/nlmodels/median/median.mod. A GAMS translation is:

$set m 19
 
set i /i1*i%m%/;
 
parameter a(i);
a(i) = uniform(0,1);
 
variable
   x 'median'
   z 'objective'
;
 
equation sumdist;
 
sumdist.. z =e= sum(i, abs(x-a(i)));
 
x.l = 0.5;
 
model median /sumdist/;
solve median minimizing z using dnlp;

Some results:

solver	x	z	Note
Minos	0.5002	4.9424	The current point cannot be improved.
Snopt	0.5002	4.9424	The current point cannot be improved.
Conopt	0.5002	4.9424	Convergence too slow
Mosek	1.864341E+12	3.542248E+13	TERMINATED BY SOLVER
Ipopt	0.5002	4.9424	Restoration Failed!
Knitro	0.5002	-0.0002	TERMINATED BY SOLVER
Baron,optcr=0,reslim=25	0.5002	4.9424	Listing file has contradictory info
OQNLP	0.5002	4.9424	NORMAL COMPLETION

Baron mentions on the log that bounds are converging to 0.494241D+01, but in the listing file it says Best possible = -1E50.

Note: This model can be reformulated as an LP:

$set m 19
 
set i /i1*i%m%/;
 
parameter a(i);
a(i) = uniform(0,1);
 
variable
   x    'median'
   e(i) 'abs deviation'
   z    'objective'
;
 
equation
  sumdist
  e1(i)
  e2(i)
;
 
sumdist.. z =e= sum(i, e(i));
e1(i)..   -e(i) =l= x-a(i);
e2(i)..             x-a(i) =l= e(i);
 
 
model median /all/;
solve median minimizing z using lp;

Recursive Descent Parser

While implementing a recursive descent parser for Excel formula's I noticed the introduction at http://en.wikipedia.org/wiki/Recursive_descent_parser.

Erratic MIQP behavior

GAMS MIQP Model:

variable z;
binary variable x;
equation e;
e.. z =e= sqr(x);
model m/e/;
solve m minimizing z using miqcp;

Results with GAMS 22.7:

1. Cplex, SBB, DICOPT,COINBONMIN,BARON: OK
2. Xpress: GAMS link does not implement MIQCP (Xpress has MIQP capability I believe)
3. MOSEK: crashes

We'll check what happens with upcoming 22.8 version.

Excel Incompatibility

Some versions of Excel allow:

  Set ws = .Sheets.Add(After:=.Sheets(Sheets.Count), _
                 Type:="worksheet", Count:=1)

This syntax is used in Microsoft examples, see e.g.: http://support.microsoft.com/kb/165238. On some other versions this may fail however with a message: "runtime error 1004" --- "worksheet.xls can not be found". More reliable is:

Set ws = .Sheets.Add(After:=.Sheets(Sheets.Count), _
             Type:=xlWorksheet, Count:=1)

Downloaded CHM files may not work

http://blogs.microsoft.co.il/blogs/tamir/archive/2007/05/03/Why-my-CHM-file-does-not-work_3F00_.aspx

Counting to one is difficult

Running

FILE OPTFILE /conopt.opt/;
  PUT OPTFILE;
  PUTCLOSE 'lsismp = t' / 'lstria = t';
VARIABLES X;
EQUATIONS E1;
E1 ..  X * 1E-12 =E= 1E-10;
MODEL M /ALL/ ;
M.OPTFILE = 1;
X.L = 1E+2; SOLVE M USING CNS;

gives the error message:

 ** Error in Square System: The number of non-fixed variables (1)
    is not equal to the number of equations (1)

Invert: An External Matrix Inversion program for GAMS

Matrix Inversion

It is possible to perform a matrix inversion using a standard LP solver. See http://www.amsterdamoptimization.com/pdf/lineq.pdf for a detailed description and examples of how this can be achieved. Sometimes it is quicker however just to call an external program to calculate the inverse of a square matrix. The program uses LAPACK's DGESV subroutine.

Download

Attach:invert.exe

Usage

The command line parameters are explained when running invert without parameters:

C:\>invert
 ----------------------------------------
 INVERT: matrix inversion
 Erwin Kalvelagen, Amsterdam Optimization
 ----------------------------------------
 Usage
    > invert gdxin i a gdxout inva
 where
    gdxin  : name of gdxfile with matrix
    i      : name of set used in matrix
    a      : name of 2 dimensional parameter inside gdxin
    gdxout : name of gdxfile for results (inverse matrix)
    inva   : name of 2 dimensional parameter inside gdxout
  
 Calculates the inverse of a non-singular matrix a(i,j) where
 i and j are aliased sets. inva will contain the inverse
 matrix.
  
  
C:\>

Example 1

The input looks like:

$ontext
 
   Finds the inverse of a matrix through an external program
 
   Erwin Kalvelagen, march 2005
 
   Reference: model gauss.gms from the model library
              http://www.gams.com/modlib/libhtml/gauss.htm
 
$offtext
 
set i  /i1*i3 /;
alias (i,j);
 
table a(i,j) 'original matrix'
      i1     i2     i3
i1    1      2      3
i2    1      3      4
i3    1      4      3
;
 
parameter inva(i,j) 'inverse of a';
 
execute_unload 'a.gdx',i,a;
execute '=invert.exe a.gdx i a b.gdx inva';
execute_load 'b.gdx',inva;
 
display a,inva;

The program has an updated GDXIO interface that both works with old 22.5 GAMS and earlier and with new GAMS 22.6 and later. The GAMS people broke the GDX interface so these versions are really incompatible. However, some of my clients were using a combination of older and newer GAMS systems, so I had to spend lots of time to write a GDX interface that both work with all these GAMS systems. Here is a log from a run with 22.6:

--- Job invert1.gms Start 02/22/08 13:29:46
GAMS Rev 149  Copyright (C) 1987-2007 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen                               G070509/0001CE-WIN
          GAMS Development Corporation                               DC4572
--- Starting compilation
--- invert1.gms(28) 3 Mb
--- Starting execution: elapsed 0:00:00.016
--- invert1.gms(25) 4 Mb
 ----------------------------------------
 INVERT: matrix inversion
 Erwin Kalvelagen, Amsterdam Optimization
 ----------------------------------------
 DLL:GDX Library   Dec 24, 2007 WIN.FR.NA 22.6 239.000.000.vis P3PC
 GDXin:a.gdx
 Input set:i
 Input parameter:a
 LAPACK DGESV, n=           3
 GDXout:b.gdx
 Output parameter:inva
 Done
--- invert1.gms(26) 4 Mb
--- GDXin=c:\tmp\New Folder\b.gdx
--- invert1.gms(28) 4 Mb
*** Status: Normal completion
--- Job invert1.gms Stop 02/22/08 13:29:46 elapsed 0:00:00.063

Example 2

$ontext
 
   Finds the inverse of a matrix through an external program
 
   Erwin Kalvelagen, march 2005
 
   Reference: model gauss.gms from the model library
              http://www.gams.com/modlib/libhtml/gauss.htm
 
$offtext
 
set i  /i1*i5/;
alias (i,j);
 
parameter a(i,j) 'Hilbert matrix';
a(i,j) = 1/(ord(i)+ord(j)-1);
 
parameter inva(i,j) 'inverse of a';
 
execute_unload 'a.gdx',i,a;
execute '=invert.exe a.gdx i a b.gdx inva';
execute_load 'b.gdx',inva;
 
display a,inva;

The output with GAMS 22.5 is:

--- Job invert2.gms Start 02/22/08 13:32:23
GAMS Rev 148  Copyright (C) 1987-2007 GAMS Development. All rights reserved
Licensee: Erwin Kalvelagen                               G070509/0001CE-WIN
          GAMS Development Corporation                               DC4572
--- Starting compilation
--- invert2.gms(24) 3 Mb
--- Starting execution
--- invert2.gms(21) 4 Mb
 ----------------------------------------
 INVERT: matrix inversion
 Erwin Kalvelagen, Amsterdam Optimization
 ----------------------------------------
 DLL:_GAMS_GDX_237_2007-01-09
 GDXin:a.gdx
 Input set:i
 Input parameter:a
 LAPACK DGESV, n=           5
 GDXout:b.gdx
 Output parameter:inva
 Done
--- invert2.gms(22) 4 Mb
--- GDXin=c:\tmp\New Folder\b.gdx
--- invert2.gms(24) 4 Mb
*** Status: Normal completion
--- Job invert2.gms Stop 02/22/08 13:32:23 elapsed 0:00:00.078

Notes

• For Lapack info see: http://www.netlib.org/lapack/
• The program can efficiently invert matrices up to (1000 x 1000). For bigger matrices, the dense inversion routine becomes slower and memory requirements become large.
• Source: Main program Attach:invert.f90
  
  

User forums, mailing lists

• Cplex: http://forums.ilog.com/optimization/index.php
• Xpress: http://groups.google.com/group/xpressmp
• AMPL: http://groups.google.com/group/ampl
• Lp_solve: http://tech.groups.yahoo.com/group/lp_solve/
• glpk: http://lists.gnu.org/archive/html/help-glpk/
• Coin or: http://list.coin-or.org/mailman/listinfo
• Sci.op-research: http://groups.google.com/group/sci.op-research/

Update:

• AIMMS: http://groups.google.com/group/aimms
• KNITRO: http://groups.google.com/group/KNITRO
• Xpress new forum: http://discuss.fico.com/blaze/board?board.id=xpressdiscuss

GAMS external program cannot write to listing file

A GAMS external program (invoked by $call or execute during respectively compile time and execution time) cannot write messages to the listing file. It would be useful to have this capability.

Delphi foreach (COM) implementation

Comlib (http://www.techvanguards.com/com/resources/downloads.asp) has an implementation of the foreach construct (used in COM Automation). Example:

var
  r : OleVariant; // range
  c : OleVariant; // cell
  Enum : TEnumVariant;

begin
   Enum := TEnumVariant.Create(r);
   while(Enum.ForEach(c)) do begin
     ....do something with c...
   end;
end;

Some minor edits are needed for Delphi 2007.

CSV File Reference Info

http://en.wikipedia.org/wiki/Comma-separated_values

Second derivative timings w. GAMS/CONOPT

The AIMMS system (a competitor of GAMS and AMPL) has according to their website much faster second derivative evaluation than GAMS: "The AIMMS implementation of the Hessian computation is extremely efficient, with reported speed ups of upto a factor 10 compared to the GAMS implementation." (From http://www.aimms.com/aimms/product/modeling_language.html). In GAMS/CONOPT the timings are reported in the listing file:

 CONOPT time Total                            0.037 seconds
   of which: Function evaluations             0.000 =  0.0%
             1st Derivative evaluations       0.000 =  0.0%
             2nd Derivative evaluations       0.000 =  0.0%
             Directional 2nd Derivative       0.000 =  0.0%

My models tend to show a relatively small percentage of time spent in the 2nd derivatives (not sure if you can conclude from this I am a good modeler). But if you have a model that shows poor timings here, I am sure the GAMS people would like to know about this, so send the model to support@gams.com. Furthermore, in such a case you may want to check if AIMMS does better on your model.

Lahey Fortran division by zero

It looks like the Lahey/Fujitsu fortran compiler does not trap divisions by zero. The code will continue with a NaN value. To generate traps add the compiler flag: -trap d.

Reading Excel Spreadsheets from Delphi

The basic code looks like:

var
  XlsObj : OleVariant;
  BookObj : OleVariant;
  SheetObj : OleVariant;
  
begin
  XlsObj := CreateOleObject('Excel.Application');
  XlsObj.DisplayAlerts := false;  // prevent warnings from popping up
  BookObj := XlsObj.Workbooks.Open(XlsFileName,,true); // readonly
  SheetObj := BookObj.Sheets[SheetName];
  do_some_work;
  BookObj.Close(false);
end;

Rolling horizon implementation in GAMS

A rolling horizon algorithm can be used to speed up finding a good integer solution. Instead of solving the model for all time periods at once we split the horizon into slices. To mitigate end-of-horizon effects, we can overlap the slices. Here is a picture of a power planning model:

(click to enlarge). Instead of solving for all periods at once we solve here in 5 steps. In each step the green part indicates binary variables are relaxed (continuous between zero and one). The yellow part indicates integer variables that are being solved for. The blue part indicates integer variables fixed to their previously found solution. Optionally we can try to improve the solution using a complete model (this is the grey part in the picture).

The GAMS code to implement this can look like:

sets
  subiter 'rolling horizon iteration' /iter1*iter5/
  relaxed(subiter,yr) /
      iter1.(2018*2037)
      iter2.(2023*2037)
      iter3.(2028*2037)
      iter4.(2033*2037)
   /
  fixed(subiter,yr) /
      iter2.(2007*2012)
      iter3.(2007*2017)
      iter4.(2007*2022)
      iter5.(2007*2027)
  /
;
  
*   Solve GEM:
    gem.optfile=1;
    gem.reslim=1000;
    gem.optcr=0;
    gem.optca=0;
    loop(subiter,
        GENBLDINT.prior(s,yr) = 1;
        loop(relaxed(subiter,yr),
           GENBLDINT.prior(s,yr) = INF;
        );
        loop(fixed(subiter,yr),
           GENBLDINT.fx(s,yr) = GENBLDINT.L(s,yr);
        );
        SOLVE GEM USING MIP MINIMIZING TOTALCOST ;
        gem.optfile=2;
    );
    gem.optfile=3;
    gem.reslim=10000;
    GENBLDINT.prior(s,yr) = 1;
    GENBLDINT.lo(s,yr) = 0;
    GENBLDINT.up(s,yr) = 1;
    SOLVE GEM USING MIP MINIMIZING TOTALCOST ;

Here we use the trick that X.PRIOR=INF means that X is relaxed. The final model uses the Cplex MIPSTART option to branch towards the solution found by the rolling horizon algorithm.

Big-M constraints in portfolio models

In general it is important to make big-M constants as small as possible.


> But, the sample i post is really a simplification for the question.
> In some cases, i must have a value for BigM greater than 100.000,
> in order to treat the real problem. For example, if [Rn] represent
> possibles ensures for diferent investiments. In this case, possible
> values for any [Rn] can be bigger than 1.000.000 ($).
>
> If we have this situation, what is a possible strategy to treat this
> question ? Maybe use a "alternative" scale for values, dividing every
> value by a common factor (p.e: /1000) ?


In portfolio models you can often express your decision variables in
terms of fractions, or in different words, use a budget of $1 instead
of $1,000,000. In the report writing code just multiply the fractions
by the actual budget.

Otherwise look into Indicator Constraints as offered by Cplex.

GAMS/IDE wish

Currently Ctrl-Z can only undo up to the last time the Run button was pressed. In some cases it can be helpful to provide undo support also for changes before that point.

In search of missing empty lines

I looked at a GAMS listing file that essentially looks like:

----     17 aaa
 
----     18 aaa
 
----     19 aaa
----     21 aaa
 
----     22 aaa
 
----     23 aaa

Why did some entries have a missing empty line between them? Not really a useful way to spent my time, but I wanted to know the answer. The reason turned out to be a solve statement with all output turned off. This does not only not write anything, it actually eats up a line from the listing file! Here is the code to demonstrate:

variable x;
equation e;
e.. x=e=1;
model m/e/;
 
*
* options to reduce output and speed up things
*
option limrow=0;
option limcol=0;
m.solprint=2;
m.solvelink=2;
 
display "aaa";
display "aaa";
display "aaa";
solve m using lp minimizing x;
display "aaa";
display "aaa";
display "aaa";

I believe this is only cosmetic and that no lines with valuable information are lost.

memory availability under GAMS

Under Windows we have the following possible cases:

1. 32 bit GAMS under 32 bit Windows gives you a 2 GB process limit in theory (in practice this is 1.6 GB).
2. 32 bit GAMS under 32 bit Windows with /3GB flag. It is possible to configure Windows XP to get you 3 GB per process (a little bit less in practice). Here is more info: http://www.microsoft.com/whdc/system/platform/server/PAE/PAEmem.mspx and http://technet.microsoft.com/en-us/library/bb124810.aspx. On Vista this is different: you need to use BCDEDIT as follows (you need administrator rights):
   bcdedit /set IncreaseUserVa 3072
   See http://www.prophotowiki.com/w/index.php/IncreaseUserVa for notes. Warning: you may want to ask your system administrator to do this for you.
3. 32 bit GAMS under 64 bit Windows. This gives you close to 4 GB.
4. 64 bit GAMS under 64 bit Windows. This gives you even more. Some NLP solvers still have a 8GB limit but otherwise you can go much higher. Note that 64 bit GAMS is using more memory than 32 bit GAMS (pointers are 64 bit instead of 32 bit). So if you only have 3 or 4 GB RAM it does not buy you much using the 64 bit GAMS version.

GAMS/Xpress bug

It is not clear why GAMS/Xpress terminated here.

               S O L V E      S U M M A R Y
 
     MODEL   lp_model            OBJECTIVE  objvar
     TYPE    MIP                 DIRECTION  MAXIMIZE
     SOLVER  XPRESS              FROM LINE  53860
 
**** SOLVER STATUS     4 TERMINATED BY SOLVER     
**** MODEL STATUS      9 INTERMEDIATE NON-INTEGER 
**** OBJECTIVE VALUE              144.0000
 
 RESOURCE USAGE, LIMIT       1612.694     10000.000
 ITERATION COUNT, LIMIT     12233      10000000
Reading parameter(s) from "C:\projects\lpsolve\xpress.opt"
>>  threads 4
>>  cutstrategy 0
>>  defaultalg 3
>>  mippresolve 011
Finished reading from "C:\projects\lpsolve\xpress.opt"
 
Xpress-MP        May  1, 2008 22.7.2 WIN 3906.4799 VIS x86/MS Windows
Xpress-MP licensed by Dash to GAMS Development Corp. for GAMS
 
The objective is a function of discrete variables only: jumpDist 1.
As a result, we set the MIP cutoff to 0.9995
Type of presolve applied: MIP presolve.
LP relaxation solved: objective =          144
Global search incomplete - no integer solution has been found
 

This may be a bug in the link instead of Xpress per se. Workaround: run again with default options will probably help.

GAMS bug: interrupt does not work

The model below does not terminate. One would expect the Interrupt button would allow you to stop the GAMS run. Unfortunately that does not work. Only the Stop button works but that stops the process very abruptly leaving an unfinished listing file.

scalar s /1/;
  
while(1,
  s = s + 1;
  s = s - 1;
);

GAMS Feature Request

For the Least Squares solvers:

• http://www.amsterdamoptimization.com/pdf/regression.pdf
• http://www.amsterdamoptimization.com/pdf/nlregression.pdf

we currently add a dummy objective to an overdetermined system of linear/non-linear equations. It would be nice if GAMS has a model type say LINSYS and NLSYS that allows to express such a system of equations. There is CNS, but that requires a square system of equations.

How to interrupt GAMS from VBA

> how can I interrupt GAMS cleanly from VBA (I call GAMS
> from VBA using CreateProcess).

You need to send a Ctrl-C to the GAMS child-process. It turns out that stuffing Ctrl-C in the keyboard buffer is unreliable. The easiest is to use a small DLL available from me that can be called as follows:

Private Declare Function GamsInterrupt Lib "gamsinterrupt.dll" _  
   (ByVal pid As Long, ByVal InterruptType As Long) As Long
'----------------------------------------------------
Private Sub SendInterrupt(InterruptType As Long)   
   Dim rc As Long   
   Dim pid As Long   
   pid = getpid()   
   rc = GamsInterrupt(ByVal pid, ByVal InterruptType)
End Sub
'----------------------------------------------------
Private Sub InterruptButton_Click()
' try to interrupt 
   Call SendInterrupt(0)
End Sub
'----------------------------------------------------
Private Sub KillButton_Click()
' try to kill 
   Call SendInterrupt(1)
End Sub

The subroutine InterruptButton_Click() mimics the Interrupt button in the GAMS/IDE and KillButton_Click() works like the Stop button. The DLL performas the low-level code-injection and remote thread managament. See also
http://www.amsterdamoptimization.com/packaging.html

A puzzling question

http://groups.google.com/group/sci.op-research/browse_thread/thread/e69aceb123a21ea5

I was a little bit too quick in responding (i.e. before I knew what the problem was). I think he wants to model

x+y=1
x,y binary

linearly but without integer variables. One can get rid of one binary variable:

y = 1-x
x binary

but if we were to be able to get rid of the integer restriction on x also, we would not need MIP solvers anymore.

A variant of the social golfer problem

http://groups.google.com/group/sci.op-research/browse_thread/thread/c45387ea26c8bf6c

> I have 28 golfers playing over four days in foursomes. I want the best
> combination. I realize not everyone can play together. Anyone with
> help? Thanks

See also the social golfer problem http://www.cs.brown.edu/~sello/golf.html.The difficult problem mentioned there has been solved using somealternative mathematical method (I forgot the details).

Analog to this problem we can minimize the number of times player A can meet player B. If this becomes 1 no player meet another player more than once.

I formulated this quickly in GAMS:

$ontext
 
 I have 28 golfers playing over four days in foursomes. I want the best
 combination. I realize not everyone can play together. Anyone with
 help? Thanks
 
 
  References:
     - problem 10 in CSPLIB (www.csplib.org)
     - Barbara Smith, "Reducing Symmetry in a Combinatorial
       Design Problem", Tech. Report, School of Computing and
       Mathematics, University of Huddersfield, 2001.
 
 
$offtext
 
 
set
   t  'days'   /day1*day4/
   i  'golfers' /golfer1*golfer28/
   g  'group'   /group1*group7/
;
 
alias(i,j);
 
binary variables    x(i,g,t)      'schedule';
positive variable   meet(i,j,g,t) 'golfer i and j meet';
free variables      maxmeet       'max number of times players meet';
 
equations
    game(i,t)          'each golver plays one game per week'
    fourplayer(g,t)    'four players per game'
    multiply1(i,j,g,t) 'linearization of multiplication (not used)'
    multiply2(i,j,g,t) 'linearization of multiplication (not used)'
    multiply3(i,j,g,t) 'linearization of multiplication'
    meetcount(i,j)
;
 
set ij(i,j);
ij(i,j)$(ord(i)>ord(j)) = yes;
 
 
*
* golfer plays one game per week
*
game(i,t).. sum(g, x(i,g,t)) =e= 1;
 
*
* four players per game
*
fourplayer(g,t).. sum(i, x(i,g,t)) =e= 4;
 
*
* linearization of x(i,g,t)*x(j,g,t)
* Note: we can relax the first two equations multiply1, multiply2
*
multiply1(ij(i,j),g,t).. meet(ij,g,t) =l= x(i,g,t);
multiply2(ij(i,j),g,t).. meet(ij,g,t) =l= x(j,g,t);
multiply3(ij(i,j),g,t).. meet(ij,g,t) =g= x(i,g,t)+x(j,g,t)-1;
 
meet.lo(ij,g,t) = 0;
meet.up(ij,g,t) = 1;
 
*
* players i and j can meet only once
*
meetcount(ij(i,j)).. sum((g,t), meet(ij,g,t)) =l= maxmeet;
 
 
 
*
* fix first round
*
set first(i,g) /
  (golfer1*golfer4).group1
  (golfer5*golfer8).group2
  (golfer9*golfer12).group3
  (golfer13*golfer16).group4
  (golfer17*golfer20).group5
  (golfer21*golfer24).group6
  (golfer25*golfer28).group7
/;
x.fx(first,'day1') = 1;
 
model m /game,fourplayer,multiply3,meetcount/;
solve m minimizing maxmeet using mip;
 
* check
parameter meetcount2(i,j) "sanity check";
meetcount2(i,j)$(not sameas(i,j)) = round(sum((g,t),x.l(i,g,t)*x.l(j,g,t)));
option meetcount2:0:1:1;
display meetcount2;
 
options x:0:2:1;
display x.l;
 

This gave a solution within 11 seconds using Cplex:

              S O L V E      S U M M A R Y
 
    MODEL   m                   OBJECTIVE  maxmeet
    TYPE    MIP                 DIRECTION  MINIMIZE
    SOLVER  CPLEX               FROM LINE  87
 
**** SOLVER STATUS     1 NORMAL COMPLETION
**** MODEL STATUS      1 OPTIMAL
**** OBJECTIVE VALUE                1.0000
 
 RESOURCE USAGE, LIMIT         10.924      1000.000
 ITERATION COUNT, LIMIT     28252         10000

with the following schedule:

----     96 VARIABLE x.L  schedule
 
                      day1        day2        day3        day4
 
golfer1 .group1           1                                   1
golfer1 .group5                       1           1
golfer2 .group1           1           1
golfer2 .group3                                   1
golfer2 .group5                                               1
golfer3 .group1           1
golfer3 .group2                                               1
golfer3 .group4                       1
golfer3 .group7                                   1
golfer4 .group1           1
golfer4 .group4                                   1
golfer4 .group6                       1
golfer4 .group7                                               1
golfer5 .group2           1
golfer5 .group5                       1
golfer5 .group7                                   1           1
golfer6 .group1                                               1
golfer6 .group2           1           1           1
golfer7 .group1                       1
golfer7 .group2           1
golfer7 .group3                                               1
golfer7 .group6                                   1
golfer8 .group2           1
golfer8 .group4                                   1           1
golfer8 .group7                       1
golfer9 .group3           1           1
golfer9 .group4                                               1
golfer9 .group7                                   1
golfer10.group1                                   1
golfer10.group2                       1
golfer10.group3           1
golfer10.group7                                               1
golfer11.group3           1                       1
golfer11.group6                                               1
golfer11.group7                       1
golfer12.group3           1                                   1
golfer12.group4                       1           1
golfer13.group4           1
golfer13.group6                       1           1           1
golfer14.group3                       1           1
golfer14.group4           1
golfer14.group7                                               1
golfer15.group1                       1           1
golfer15.group4           1                                   1
golfer16.group2                       1
golfer16.group4           1
golfer16.group5                                   1           1
golfer17.group3                                   1
golfer17.group4                                               1
golfer17.group5           1           1
golfer18.group1                                   1
golfer18.group3                                               1
golfer18.group5           1
golfer18.group6                       1
golfer19.group2                                   1
golfer19.group4                       1
golfer19.group5           1
golfer19.group6                                               1
golfer20.group2                                               1
golfer20.group5           1                       1
golfer20.group7                       1
golfer21.group2                                   1
golfer21.group3                       1                       1
golfer21.group6           1
golfer22.group2                       1                       1
golfer22.group6           1                       1
golfer23.group1                                   1           1
golfer23.group4                       1
golfer23.group6           1
golfer24.group4                                   1
golfer24.group5                       1                       1
golfer24.group6           1
golfer25.group3                       1
golfer25.group5                                               1
golfer25.group6                                   1
golfer25.group7           1
golfer26.group1                       1
golfer26.group5                                   1
golfer26.group6                                               1
golfer26.group7           1
golfer27.group1                                               1
golfer27.group7           1           1           1
golfer28.group2                                   1           1
golfer28.group6                       1
golfer28.group7           1

GAMS/Cplex docs headers are messed up

The GAMS/Cplex docs contain some section headers that are messed up. (Actually with LaTeX this usually works correctly; you really have to mess up the input to achieve this).

I don't think any important information is missing.

Modeling question

http://tech.groups.yahoo.com/group/lp_solve/message/10534.

Looks to me using some homotopy method for this is inappropriate (I thought these were for nonlinear problems instead of discrete problems).

GAMS/IDE suggestion

It would be helpful if the edit boxes in the Find/Replace window would be larger if the window is enlarged. Long paths are not completely visible using the default sizes for the edit boxes. Click on the picture to see what I mean. The edit boxes should/could scale with the window size.

GAMS/IDE bug

Ctrl-T eats up too many characters:

›  // hello
› -10

In both cases Ctrl-T should remove blanks up-to the first non-blank character. However it will also remove the // and - characters.

If this happens, use Ctrl-Z to undo.

MWsnap for making screen shots

http://www.mirekw.com/winfreeware/mwsnap.html is a useful tool to create screen shots. I use it a lot.

GAMS/IDE bug

F3 does not repeat a search/replace operation (it becomes a search only). Sometimes it is useful to verify each search/replace and F3 can be used to do this. However, in the GAMS/IDE F3 does not repeat a complete search/replace operation.

Demonstration of numerical issues in LS calculations.

Probably the worst algorithm to perform a least squares estimation is to evaluate

b = (X'X)^-1X'y

This is illustrated in http://amsterdamoptimization.com/models/regression/longley2.gms. The solution compares as follows:

 parameter  reference                 LS solver     b=inv(X'X)X'y
    B0     -3482258.63459582         -3.48226E+6       65317.000
    B1      15.0618722713733         15.06187227          0
    B2     -0.358191792925910E-01    -0.03581918          0
    B3     -2.02022980381683         -2.02022980          0
    B4     -1.03322686717359         -1.03322687          0
    B5     -0.511041056535807E-01    -0.05110411          0
    B6      1829.15146461355         1.829151E+3          0

A numerically sound method is documented in http://amsterdamoptimization.com/pdf/regression.pdf. The above model can then be coded as: http://amsterdamoptimization.com/models/regression/longley.gms. The results are listed in the column "LS Solver".

Excel PivotTable using external data

When storing a GAMS symbol each nonzero in a separate row in Excel, we quickly run out of space on older versions of Excel. Excel before Office 2007 could only handle up to 65,536 rows. To display larger pivot tables, we can write the GAMS symbol to a CSV file and use this CSV file as external data store for the Excel pivot table. This allows for larger data sets than the above mentioned row limit. See also http://www.amsterdamoptimization.com/pdf/toolhelp.pdf.

AMPL defined variables

AMPL has a construct called defined variables that can sometimes reduce the size of an NLP by a large amount. Here are some statistics for a client model (confidential so I cannot share it):

GAMS:

BLOCKS OF EQUATIONS           6     SINGLE EQUATIONS       26,192
BLOCKS OF VARIABLES           5     SINGLE VARIABLES       26,201
NON ZERO ELEMENTS       119,031     NON LINEAR N-Z         45,220
DERIVATIVE POOL               8     CONSTANT POOL           1,541
CODE LENGTH             359,381

AMPL:

Substitution eliminates 26180 variables.
Adjusted problem:
10 variables, all nonlinear
1 constraint, all linear; 10 nonzeros 
1 nonlinear objective; 10 nonzeros.

GAMS has some similar facilities built in the LGO solver, but these are not available as a general feature.

Loops over different data sets

See http://groups.google.com/group/sci.op-research/browse_thread/thread/78456e35ddfc7141. This is actually a question that is often asked. The answer is really dependent on the situation. In this case the single case digests the input during compile time (using $include). This can not be coded directly by placing a loop around it (a loop is runtime). In this case it would probably the easiest the implement the single case as a complete GAMS file passing on the name of the input file and the output file as command line parameters. I.e.

* %input% is input file name (XLS file)
* %output% is output gdx file

$if not set input $set input E:\simuldata1\datag1.xls 
$if not set output $set output E:\simuldata\data1.xls  
 
Set s1  /
$call =xls2gms "i=%input%" o=E:\work\dataif1.inc R=datag1!a2:a12
$include E:\work\dataif1.inc
/;
  
Parameter p(s1) /
$call =xls2gms "i=%input%" o=E:\work\data21.inc R=datag1!j2:k12
$include E:\work\data2l.inc
/;
 
......
 
Execute_Unload "data2.gdx" V1
Execute 'GDXXRW.EXE data2.gdx sq=0 o=%output% var=V1 rng=sheet2!e2 rdim=1';

Note: the $if statements allow the model to be run outside the loop as well.

Now generate the loop in a separate GAMS file:

set i /1*100/;
 
file f /x.cmd/;
 
loop(i,
   put f,"gams singlecase "
         "--input=E:\simuldata" i.tl:0 "\datag" i.tl:0 ".xls "
         "--output=E:\simuldata\data" i.tl:0 ".xls "
         "O=E:\simuldata\singlecase" i.tl:0 ".lst"/;
);
putclose f;
  
execute "x.cmd";

This will generate a batch file:

gams singlecase --input=E:\simuldata1\datag1.xls --output=E:\simuldata\data1.xls O=E:\simuldata\singlecase1.lst
gams singlecase --input=E:\simuldata2\datag2.xls --output=E:\simuldata\data2.xls O=E:\simuldata\singlecase2.lst
gams singlecase --input=E:\simuldata3\datag3.xls --output=E:\simuldata\data3.xls O=E:\simuldata\singlecase3.lst
gams singlecase --input=E:\simuldata4\datag4.xls --output=E:\simuldata\data4.xls O=E:\simuldata\singlecase4.lst
gams singlecase --input=E:\simuldata5\datag5.xls --output=E:\simuldata\data5.xls O=E:\simuldata\singlecase5.lst
....

Instead of writing a batch file it is also possible to directly generate a command line using the infamous put_utility command:

put_utility "shell" / "gams singlecase --input=E:\simuldata" i.tl:0 "\datag" i.tl:0 ".xls --output=E:\simuldata\data" i.tl:0 ".xls O=E:\simuldata\singlecase" i.tl:0 ".lst";

Probably the best is to consult support@gams.com as this is very hairy.

Often it is better to reorganize the data flow and read all data in the beginning, then perform a loop to calculate the results, and then export all results. The input data and results have then an extra index.

Any uglier syntax possible?

To compose a file name for a gdx file, one can do:

set mds /mds1*mds3/;
parameter x;
 
* needed for put_utility
file dummy;
put dummy;
 
* remove all gdx files
execute 'del /q *.gdx';
 
loop(mds,
 
* solve etc here

* say solution is x
    x = ord(mds);
    display x;
 
* filename := 'sol_mds1',  'sol_mds2', ....
    put_utility 'gdxout' / 'sol_' mds.tl:0;
 
* write to gdx file with changed name
    execute_unload x;
);
 
* show what gdx files we now have
execute 'dir *.gdx';
 
* merge gdx files
execute 'gdxmerge sol_*.gdx';
 
* show that it created merged.gdx
execute 'dir *.gdx';
 
* show content of merged.gdx
execute 'gdxdump merged.gdx';

In my experience this is very difficult to explain. The put_utility syntax should really be replaced by something more elegant (not a very high bar).

Notice that we need a put file around even if we don't use it. It should not be too difficult to hide this from the modeler and handle this inside GAMS.

Gamma Distribution CDF

See http://groups.google.com/group/ampl/browse_thread/thread/10861bd79114dfec. In GAMS we have the GAMMAREG function (see http://www.gams.com/~erwin/specfun.pdf) so no need to program this yourself using external functions.

The Gamma CDF from http://en.wikipedia.org/wiki/Gamma_distribution denoted by F(x;k,θ) with

• scale parameter θ
• shape parameter k
• mean: k·θ
• variance: k·θ²

can be expressed in GAMS by:

F(x;k,θ) = GAMMAREG(x/θ,k)

Both first and second derivatives are available for solvers calling this function.

More execute_unload

This seems allowed but has unexpected results.

scalar a /1/; 
execute_unload a;

This will go to a gdx file 'xllink.gdx'???? Be warned: this construct will overwrite any file xllink.gdx even though you did not specify this filename in the GAMS model (luckily chance is pretty small).

GAMS: execute_unload

It would be useful to be able to write to the same GDX file in different steps. I.e. something like:

execute_unload "a.gdx",a;

a(i) = a(i)/1000;
b(i) = c(i) + a(i);

execute_unload "a.gdx","append",b; 

Currently we cannot do this: there is no facility to append.

GAMS/IDE suggestion

The LST viewer can be improved by allowing easy navigation to syntax errors in the listing file. In many cases you don't even need to consult the listing file to fix errors using the IDE, but sometimes it is needed to inspect the listing file. (The previous is an example). This may help to find syntax error messages in the listing file (especially when there are comments with stars so searching for **** is not convenient).