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 ;

When do two intervals overlap?

For a scheduling tool I encountered the following problem. For a set of intervals generate a list of a pairs that overlap. For those pairs (A,B) I can generate a constraint (Note: we are smart and do not generate both x_A+x_B≤1 and x_B+x_A≤1) The test if two intervals A=[a₁,a₂], B=[b₁,b₂] overlap is:

or

depending on borderline cases (the last version is called the “strict” version in the following). Here is the demonstration:

Note that if the intervals are determined endogenously in the model, we often write “no-overlap” constraints for two intervals A=[a₁,a₂], B=[b₁,b₂] as:

I.e. B starts after A ends or A starts after B ends. Indeed, we can see that

which derives our second “strict” form of the overlap conditions.

Coin CBC command line

For a small project I was trying out the binary distribution of COIN-OR’s CBC. The good news, it can solve the problem quickly (provided I allow a  reasonable gap). I encountered a few issues:

1. The Windows executables are not built with parallel threading support.
2. The CSV writer is too simplistic. Names (with commas in them) are not properly quoted.

Threads option ignored

C:\projects\test>cbc model1a.lp ratio 0.05 threads 8 solve printing csv solution model1a.csv Welcome to the CBC MILP Solver Version: 2.8.8 Build Date: Dec 23 2013 Revision Number: 1998 command line - cbc model1a.lp ratio 0.05 threads 8 solve printing csv solution model1a.csv (default strategy 1) ratioGap was changed from 0 to 0.05 No match for threads - ? for list of commands No match for 8 - ? for list of commands Continuous objective value is 57.0769 - 0.02 seconds Cgl0002I 30 variables fixed Cgl0003I 0 fixed, 0 tightened bounds, 527 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 475 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 429 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 291 strengthened rows, 0 substitutions Cgl0003I 0 fixed, 0 tightened bounds, 69 strengthened rows, 0 substitutions Cgl0004I processed model has 743 rows, 499 columns (346 integer) and 2429 elements Cbc0038I Pass   1: suminf.    0.00000 (0) obj. 12146 iterations 431 Cbc0038I Solution found of 12146 Cbc0038I Before mini branch and bound, 316 integers at bound fixed and 140 continuous Cbc0038I Full problem 743 rows 499 columns, reduced to 247 rows 145 columns Cbc0038I Mini branch and bound improved solution from 12146 to 166 (0.50 seconds) 0  Obj 166 Dual inf 0.22137435 (11) End of values pass after 26 iterations 26  Obj 166 Dual inf 549.2981 (110) Perturbing problem by 0.001 % of 0.92934395 - largest nonzero change 1.9181672e-005 (% 0.0011298852) - largest zero chan ge 1.9246056e-005 120  Obj 166 Optimal - objective value 166 Cbc0038I Round again with cutoff of 164.5 Cbc0038I Pass   2: suminf.    7.33333 (22) obj. 164.5 iterations 163 Cbc0038I Pass   3: suminf.    6.22857 (20) obj. 164.5 iterations 71 Cbc0038I Pass   4: suminf.    5.33336 (18) obj. 164.5 iterations 34 Cbc0038I Pass   5: suminf.    4.10000 (14) obj. 164.5 iterations 28 Cbc0038I Pass   6: suminf.    0.00750 (3) obj. 164.5 iterations 166 Cbc0038I Solution found of 162 Cbc0038I Before mini branch and bound, 311 integers at bound fixed and 141 continuous Cbc0038I Full problem 743 rows 499 columns, reduced to 296 rows 158 columns Cbc0038I Mini branch and bound did not improve solution (0.74 seconds) Cbc0038I Round again with cutoff of 159.8 Cbc0038I Pass   7: suminf.    8.02667 (26) obj. 159.8 iterations 29 Cbc0038I Pass   8: suminf.    7.10546 (36) obj. 159.8 iterations 122 Cbc0038I Pass   9: suminf.    7.12365 (46) obj. 159.8 iterations 40 Cbc0038I Pass  10: suminf.   18.70686 (90) obj. 159.8 iterations 281 Cbc0038I Pass  11: suminf.    9.66191 (71) obj. 159.8 iterations 220 Cbc0038I Pass  12: suminf.    9.07849 (73) obj. 159.8 iterations 59 Cbc0038I Pass  13: suminf.    9.67900 (63) obj. 159.8 iterations 157 Cbc0038I Pass  14: suminf.    8.05376 (64) obj. 159.8 iterations 92 Cbc0038I Pass  15: suminf.    9.90251 (61) obj. 159.8 iterations 100 Cbc0038I Pass  16: suminf.    8.87667 (61) obj. 159.8 iterations 59 Cbc0038I Pass  17: suminf.   11.47273 (46) obj. 159.8 iterations 120 Cbc0038I Pass  18: suminf.   10.18668 (44) obj. 159.8 iterations 33 Cbc0038I Pass  19: suminf.    7.58501 (60) obj. 159.8 iterations 60 Cbc0038I Pass  20: suminf.    6.99899 (78) obj. 159.8 iterations 36 Cbc0038I Pass  21: suminf.    7.66668 (32) obj. 159.8 iterations 128 Cbc0038I Pass  22: suminf.    7.13455 (64) obj. 159.8 iterations 106 Cbc0038I Pass  23: suminf.    7.20358 (59) obj. 159.8 iterations 80 Cbc0038I Pass  24: suminf.    7.13455 (64) obj. 159.8 iterations 67 Cbc0038I Pass  25: suminf.    8.20892 (48) obj. 159.8 iterations 257 Cbc0038I Pass  26: suminf.    7.32000 (30) obj. 159.8 iterations 124 Cbc0038I Pass  27: suminf.    9.03479 (27) obj. 159.8 iterations 92 Cbc0038I Pass  28: suminf.    7.34334 (33) obj. 159.8 iterations 59 Cbc0038I Pass  29: suminf.    7.32000 (30) obj. 159.8 iterations 53 Cbc0038I Pass  30: suminf.    7.32000 (30) obj. 159.8 iterations 31 Cbc0038I Pass  31: suminf.   16.85062 (65) obj. 159.8 iterations 288 Cbc0038I Pass  32: suminf.    5.71250 (21) obj. 159.8 iterations 277 Cbc0038I Pass  33: suminf.    5.00000 (15) obj. 159.8 iterations 69 Cbc0038I Pass  34: suminf.    5.00000 (15) obj. 159.8 iterations 72 Cbc0038I Pass  35: suminf.    5.00000 (15) obj. 159.8 iterations 136 Cbc0038I Pass  36: suminf.    9.86669 (43) obj. 159.8 iterations 256 Cbc0038I No solution found this major pass Cbc0038I Before mini branch and bound, 178 integers at bound fixed and 89 continuous Cbc0038I Full problem 743 rows 499 columns, reduced to 542 rows 307 columns - 15 fixed gives 331, 197 - ok now Cbc0038I Full problem 743 rows 499 columns, reduced to 328 rows 197 columns Cbc0038I Mini branch and bound did not improve solution (1.22 seconds) Cbc0038I After 1.23 seconds - Feasibility pump exiting with objective of 162 - took 0.92 seconds Cbc0012I Integer solution of 162 found by feasibility pump after 0 iterations and 0 nodes (1.24 seconds) Cbc0038I Full problem 743 rows 499 columns, reduced to 317 rows 154 columns Cbc0038I Full problem 743 rows 499 columns, reduced to 0 rows 0 columns Cbc0031I 19 added rows had average density of 171.78947 Cbc0013I At root node, 19 cuts changed objective from 151 to 154.55353 in 100 passes Cbc0014I Cut generator 0 (Probing) - 207 row cuts average 2.5 elements, 7 column cuts (7 active)  in 0.203 seconds - new frequency is 1 Cbc0014I Cut generator 1 (Gomory) - 1562 row cuts average 349.6 elements, 0 column cuts (0 active)  in 0.842 seconds - n ew frequency is 1 Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.094 seconds - new frequency is -100 Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.037 seconds - new fr equency is -100 Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 1.313 s econds - new frequency is -100 Cbc0014I Cut generator 6 (TwoMirCuts) - 191 row cuts average 191.0 elements, 0 column cuts (0 active)  in 0.253 seconds - new frequency is 1 Cbc0011I Exiting as integer gap of 7.4464682 less than 1e-010 or 5% Cbc0001I Search completed - best objective 162, took 17263 iterations and 0 nodes (6.98 seconds) Cbc0035I Maximum depth 0, 42 variables fixed on reduced cost Cuts at root node changed objective from 151 to 154.554 Probing was tried 100 times and created 214 cuts of which 0 were active after adding rounds of cuts (0.203 seconds) Gomory was tried 100 times and created 1562 cuts of which 0 were active after adding rounds of cuts (0.842 seconds) Knapsack was tried 100 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.094 seconds) Clique was tried 100 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.037 seconds) MixedIntegerRounding2 was tried 100 times and created 0 cuts of which 0 were active after adding rounds of cuts (1.313 s econds) FlowCover was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 100 times and created 191 cuts of which 0 were active after adding rounds of cuts (0.253 seconds) Result - Optimal solution found (within gap tolerance) Objective value:                162.00000000 Lower bound:                    154.554 Gap:                            0.05 Enumerated nodes:               0 Total iterations:               17263 Time (CPU seconds):             7.05 Time (Wallclock seconds):       7.05 Option for printingOptions changed from normal to csv Total time (CPU seconds):       7.08   (Wallclock seconds):       7.08 C:\projects\test>

A better error message could indicate this executable was not compiled with parallel threading support.
I think parallel CBC does not work very well on Windows. I have tried also a GAMS version of the model using GAMS/CBC with threads=8. This leads to a crash:



Looks like this is just not working very well. May be it is good threads are disabled in the Windows executable from the binary dsitribution.

Incorrect CSV output

The column name has some embedded commas, so the name should be quoted by the CSV writer. Here is a fragment of the output.

Assign(Student1,MON,m3),0 Assign(Student1,SCH,m1),1 Assign(Student1,MON,m4),0 Assign(Student1,MON,m7),0 Assign(Student1,MON,m8),0 Assign(Student1,GES,m11),0 Assign(Student1,GES,m12),0 Assign(Student1,GES,m13),0

Update: running on Mac (GAMS/CBC)

This also looks very scary (memory overruns?):

Maximum depth 0, 6 variables fixed on reduced cost 7 added rows had average density of 51.857143 At root node, 17 cuts changed objective from 158.29268 to 161.10418 in 8 passes Cut generator 0 (Gomory) - 37 row cuts average 58.3 elements, 0 column cuts (12 active) CutS egaerncehr actoomrp l0e t(eGdo m-o rbye)s t-  o3b7j ercotwi vceu t1se +a5v0e,r atgoeo k5 82.039  eilteemreanttiso,n s0  acnodl u0m nn ocduetss  ((1182. 5a5c tsievceo nds) Maximum depth 0, 6 variables fixed on reduced cost

The CBC log is pretty dense and messy (matter of taste of course), so not always easy to see strange messages, but this gibberish is of course a sign of a problem.

Update2: Mac again

After recompiling from sources I could run the problem with multiple threads (when reading from MPS files). May be a build problem in the GAMS version of CBC. But now I see my LP files are no longer accepted:

Coin3007W ### CoinLpIO::is_invalid_name(): Name is empty Coin3007W ### CoinLpIO::are_invalid_names(): Invalid name: vnames[1396]: (null) Coin3007W ### CoinLpIO::setLpDataRowAndColNames(): Invalid row names

Not a very clear message. I would like to see what I should do to fix my file. The solution is bogus so the lp file is read in incorrectly. (To be complete: this was tar file version 2.8.12).  I am not a happy camper…