GAMS is more geared towards solving large-scale nonlinear programming problems, where just one, hopefully optimal, solution is reported. The typical NLP solvers under GAMS support this single solution paradigm (they are using numerical optimization algorithms). A symbolic math package such as Mathematica would be ideally suited for a problem of finding all roots of a polynomial.

I am a full-time consultant and provide services related to the design, implementation and deployment of mathematical programming, optimization and data-science applications. I also teach courses and workshops. Usually I cannot blog about projects I am doing, but there are many technical notes I'd like to share. Not in the least so I have an easy way to search and find them again myself. You can reach me at erwin@amsterdamoptimization.com.

## Thursday, February 26, 2009

### Finding all roots with GAMS

GAMS is more geared towards solving large-scale nonlinear programming problems, where just one, hopefully optimal, solution is reported. The typical NLP solvers under GAMS support this single solution paradigm (they are using numerical optimization algorithms). A symbolic math package such as Mathematica would be ideally suited for a problem of finding all roots of a polynomial.

## Tuesday, February 24, 2009

### Timing of Sparse Matrix Multiplication in Mathprog, Ampl and GAMS

*In Mathprog, can I use sparse matrix multplication?*

`param N := 250;`

set I := {1..N};

param A{i in I, j in I} := if (i=j) then 1;

param B{i in I, j in I} := if (i=j) then 1;

param C{i in I, j in I} := sum{k in I} A[i,k]*B[k,j];

param s := sum{i in I, j in I} C[i,j];

display s;

end;

Using the

**timeit.exe**tool from the windows 2003 resource kit, we see:

`C:\projects\glpk\glpk\bin>timeit glpsol --math x.mod`

Reading model section from x.mod...

9 lines were read

Display statement at line 7

s = 250

Model has been successfully generated

Problem has no rows

Version Number: Windows NT 6.0 (Build 6001)

Exit Time: 4:28 pm, Tuesday, February 24 2009

Elapsed Time: 0:00:51.025

Process Time: 0:00:50.076

System Calls: 3306455

Context Switches: 841545

Page Faults: 191189

Bytes Read: 36549676

Bytes Written: 1409054

Bytes Other: 120151

C:\projects\glpk\glpk\bin>timeit \projects\ampl\amplcml\ampl.exe x.mod

s = 250

Version Number: Windows NT 6.0 (Build 6001)

Exit Time: 4:28 pm, Tuesday, February 24 2009

Elapsed Time: 0:00:02.571

Process Time: 0:00:02.542

System Calls: 104208

Context Switches: 20998

Page Faults: 9767

Bytes Read: 59854

Bytes Written: 784

Bytes Other: 6240

C:\projects\glpk\glpk\bin>

I.e. Mathprog takes 50 seconds and AMPL only 3 seconds. This difference is largely due to smarter implementation in AMPL, but likely not to better sparse execution. If we fill the matrices A and B with all ones we get similar timings (50 vs 3 seconds).

When we translate this to GAMS:

`set i /1*250/;`

alias (i,j,k);

parameter A(i,j),B(i,j),C(i,j);

A(i,i) = 1;

B(i,i) = 1;

C(i,j) = sum(k, A(i,k)*B(k,j));

scalar s;

s = sum((i,j),C(i,j));

display s;

we get even better performance:

`C:\projects\glpk\glpk\bin>timeit gams.exe x.gms`

--- Job x.gms Start 02/24/09 16:33:00 WEX-WEI 23.0.2 x86_64/MS Windows

GAMS Rev 230 Copyright (C) 1987-2009 GAMS Development. All rights reserved

Licensee: Erwin Kalvelagen G090213/0001CV-WIN

Amsterdam Optimization Modeling Group DC4572

--- Starting compilation

--- x.gms(9) 3 Mb

--- Starting execution: elapsed 0:00:00.008

--- x.gms(9) 4 Mb

*** Status: Normal completion

--- Job x.gms Stop 02/24/09 16:33:00 elapsed 0:00:00.031

Version Number: Windows NT 6.0 (Build 6001)

Exit Time: 4:33 pm, Tuesday, February 24 2009

Elapsed Time: 0:00:00.109

Process Time: 0:00:00.015

System Calls: 7493

Context Switches: 1558

Page Faults: 3959

Bytes Read: 495203

Bytes Written: 418566

Bytes Other: 702554

C:\projects\glpk\glpk\bin>

I.e. only 0.1 seconds. This is indeed due to sparse processing of the matrix multiplication. If we fill A and B with ones then the GAMS timing becomes close to AMPL: about 3 seconds. Indeed when we look at the GAMS profiling information we see how sparsity is exploited in the different statements:

`**** SPARSE MODEL (Unity matrix)`

---- 1 ExecInit 0.000 0.000 SECS 3 Mb

---- 4 Assignment A 0.000 0.000 SECS 4 Mb 250

---- 5 Assignment B 0.000 0.000 SECS 4 Mb 250

---- 6 Assignment C 0.016 0.016 SECS 4 Mb 250

---- 8 Assignment s 0.000 0.016 SECS 4 Mb

---- 9 PARAMETER s = 250.000

---- 9 Display 0.000 0.016 SECS 4 Mb

**** DENSE MODEL (Matrix with ones)

---- 1 ExecInit 0.000 0.000 SECS 3 Mb

---- 4 Assignment A 0.000 0.000 SECS 6 Mb 62500

---- 5 Assignment B 0.015 0.015 SECS 8 Mb 62500

---- 6 Assignment C 2.746 2.761 SECS 13 Mb 62500

---- 8 Assignment s 0.000 2.761 SECS 13 Mb

---- 9 PARAMETER s = 1.562500E+7

---- 9 Display 0.000 2.761 SECS 13 Mb

All the action is in the assignment to C and this is where GAMS really shows it strength. To summarize, the timings are roughly:

Model | Mathprog | AMPL | GAMS |
---|---|---|---|

sparse (identity) | 50 sec | 3 sec | 0.1 sec |

dense (all ones) | 50 sec | 3 sec | 3 sec |

I suspect AMPL is not fully exploiting sparsity here. In general AMPL does however do a really good job in using sparse patterns during model generation with performance similar to GAMS. I have seen before models with Mathprog is just collapsing where the same model executes and generates very quickly under AMPL or GAMS. Nevertheless it is quite amazing how many models generate fast with Mathprog. This is likely because in many cases, although the final LP matrix is very large and sparse, substructures in the model can often be dense and relatively small. I.e. many models can be generated quite efficiently using dense technology for most operations. Only when the substructures become very large and sparse, better sparse technology -- i.e. loop over the nonzero elements only -- as used in GAMS (and AMPL) becomes necessary.

## Friday, February 20, 2009

## Wednesday, February 18, 2009

### MS Solver Foundation /Gurobi (2)

Error:

The solver(s) threw an exception while solving the model - Gurobi throttle exceeded.

Visit www.gurobi.com/html/products.html to purchase an unlock license (HostID: 8AA3C6D2).

The obj seems to be counted as an equation:

`Model[`

// generate simple MIP model with 500 variables and 500 constraints

Parameters[Integers,N=500],

Decisions[Integers[0,1],Foreach[{i,N},x[i]]],

Constraints[ Foreach[{i,N},x[i] <= 1] ],

Goals[

Maximize[Sum[{i,N},x[i]]]

//Maximize[obj->Sum[{i,N},x[i]]] use to prevent funny name for goal

]

]

does trigger the above message. From the Excel plug-in I don't see an option to use the MS MIP solver. I hope I am wrong, but I suspect there is a class of MIP models that could be solved with the previous version of the Standard Edition that can no longer be solved with the current 1.1 release.

The workaround seems obvious: purchase the enterprise edition or purchase Gurobi.

Note that the limit is on any variable, i.e. 1 integer variable + 500 continuous variables makes the model too large.

Update: If your problem is in between these limits this can help. There is a way to point the excel plugin back to the MS MIP solver through a configuration file. For more info please contact the MSF people at http://code.msdn.microsoft.com/solverfoundation. They are very helpful and responsive.

### Basis information in GAMS

` LOWER LEVEL UPPER MARGINAL`

---- EQU cost . . . 1.0000 Non-basic

---- EQU supply observe supply limit at plant i

LOWER LEVEL UPPER MARGINAL

seattle -INF 350.0000 350.0000 EPS Non-basic

san-diego -INF 550.0000 600.0000 . Basic

---- EQU demand satisfy demand at market j

LOWER LEVEL UPPER MARGINAL

new-york 325.0000 325.0000 +INF 0.2250 Non-basic

chicago 300.0000 300.0000 +INF 0.1530 Non-basic

topeka 275.0000 275.0000 +INF 0.1260 Non-basic

---- VAR x shipment quantities in cases

LOWER LEVEL UPPER MARGINAL

seattle .new-york . 50.0000 +INF . Basic

seattle .chicago . 300.0000 +INF . Basic

seattle .topeka . . +INF 0.0360 Non-basic

san-diego.new-york . 275.0000 +INF . Basic

san-diego.chicago . . +INF 0.0090 Non-basic

san-diego.topeka . 275.0000 +INF . Basic

LOWER LEVEL UPPER MARGINAL

---- VAR z -INF 153.6750 +INF . Basic

This model has 6 equations and 7 variables. So we expect 6 basic rows/columns and 7 non-basic ones. This indeed is the case.

Note: this model is

*dual degenerate*as it has an EPS marginal. This means the dual of equation supply('seattle') is numerically zero but non-basic nevertheless. This indicates presence of alternative optimal solutions.

## Friday, February 13, 2009

### MS Solver Foundation / Gurobi

## Thursday, February 12, 2009

### TSP Powerset Formulation

*> I'm a PhD student in Portugal and I´m modeling a Vehicle Routing*

> Problem. I have some difficulties to modelling this subtour

> elimination constraint (n = number of nodes):

> Problem. I have some difficulties to modelling this subtour

> elimination constraint (n = number of nodes):

*> sum(i in S, sum(j in S, x(i,j))) <= |S|-1*

> for every nonempty subset S of {2,3,..,n}

> for every nonempty subset S of {2,3,..,n}

This is the Ford, Fulkerson, Johnson (1954) formulation, based on the power of {2,3,..,n}. This generates a lot of constraints (2^(n-1)), with many not active in the optimal solution. So this approach only works for small problems.

$ontext Erwin Kalvelagen, Amsterdam OptimizationReferences:Gerhard Reinelt, Universitaet Heidelberg, TSPLIB95A.J. Orman & H.P. Williams, A Survey of Different Integer ProgrammingFormulations of the Travelling Salesman Problem$offtext $eolcom // setsi /city1*city14/ xy /x,y/ ; alias(i,j);table coordinates(i,xy) 'coordinates from tsplib'x y city1 16.47 96.10 city2 16.47 94.44 city3 20.09 92.54 city4 22.39 93.37 city5 25.23 97.24 city6 22.00 96.05 city7 20.47 97.02 city8 17.20 96.29 city9 16.30 97.38 city10 14.05 98.12 city11 16.53 97.38 city12 21.52 95.59 city13 19.41 97.13 city14 20.09 94.55 ; ** form the distance matrix*parameter latitude(i), longitude(i), degree(i,xy), minute(i,xy);degree(i,xy) = trunc(coordinates(i,xy)); minute(i,xy) = coordinates(i,xy) - trunc(coordinates(i,xy)); latitude(i) = pi*(degree(i,'x')+5.0*minute(i,'x')/3.0)/180.0; longitude(i) = pi*(degree(i,'y')+5.0*minute(i,'y')/3.0)/180.0; set nd(i,j);nd(i,j)$( not sameas(i,j)) = yes;scalar rrr /6378.388/;parameter q1(i,j),q2(i,j),q3(i,j);q1(nd(i,j)) = cos(longitude(i) - longitude(j)); q2(nd(i,j)) = cos(latitude(i) - latitude(j)); q3(nd(i,j)) = cos(latitude(i) + latitude(j)); parameter c(i,j) 'distance matrix (KM)';c(nd(i,j)) = trunc(RRR * arccos( 0.5*((1.0+q1(i,j))*q2(i,j) - (1.0-q1(i,j))*q3(i,j)) ) + 1.0); display c;** form the power set*set v(i) 'exclude city 1' /city2*city14/;abort$(card(i)<>card(v)+1) 'Set v is incorrect',i,v;set n 'subset id (max number of subsets)' /n1*n8192/;scalar nsize 'size needed for set n';nsize = 2** card(v);abort$(card(n)<nsize) 'set n is too small',nsize;setsbase 'used in next set' /no,yes/ps0(n,v,base) / system.powersetRight/ ps(n,v) 'power set'; ps(n,v) = ps0(n,v,'yes'); parameter pscount(n);pscount(n) = sum(ps(n,v),1);** form subsets, exclude subsets with cardinality = 1*set nsub(n);nsub(n)$(pscount(n)>=2) = yes;** form ps2(n,i,j) = ps(n,i) and ps(n,j)* this will simplify the subtour elimination constraint*alias(v,vv);set ps2(n,i,j);ps2(nsub(n),nd(v,vv))$(ps(n,v) and ps(n,vv)) = yes;** here is the model*binary variables x(i,j);variable z 'objective variable';equationsobj 'objective'assign1(i) 'assignment problem constraint'assign2(j) 'assignment problem constraint'subt(n) 'subtour elimination'; obj.. z =e= sum(nd(i,j), c(i,j)*x(i,j));assign1(i).. sum(nd(i,j), x(i,j)) =e= 1;assign2(j).. sum(nd(i,j), x(i,j)) =e= 1;subt(nsub(n)).. sum(ps2(n,i,j),x(i,j)) =L= pscount(n)-1;option optcr=0;model tsp /all/;solve tsp minimizing z using mip; |

Note: the construct **system.powerSetRight **is a new GAMS feature (24.0, 2013).

This solves quite fast:

--- Starting compilation IBM ILOG CPLEX 24.7.3 r58181 Released Jul 11, 2016 WEI x86 64bit/MS Windows Reading data... Iteration log . . . Nodes Cuts/ * 0 0 integral 0 3323.0000 3323.0000 45 0.00% Root node processing (before b&c): Proven optimal solution. MIP Solution: 3323.000000 (45 iterations, 0 nodes) Best possible: 3323.000000 --- Restarting execution |

##### References

- burma14.gms: above model
- burma14-3.gms: generation of powerset written in GAMS (not using
**system.PowerSetRight**). - burma14-2.gms: formulation from Svestka (J.A. Svestka,
*A continuous variable representation of the TSP*, Math Prog, v15, 1978, pp211-213) - burma14-4.gms: a cutting plane technique

### Easy MIP

`--- 122,284 rows 352,038 columns 667,797 non-zeroes`

--- 231,136 discrete-columns

Nodes Cuts/

Node Left Objective IInf Best Integer Best Node ItCnt Gap

0 0 588575.0000 606 588575.0000 14792

0 0 588575.0000 689 Cuts: 313 16547

0 0 588575.0000 736 ZeroHalf: 14 17624

* 0+ 0 587775.0000 588575.0000 17624 0.14%

0 2 588575.0000 349 587775.0000 588575.0000 17624 0.14%

Elapsed real time = 30.14 sec. (tree size = 0.00 MB, solutions = 1)

* 10+ 1 588575.0000 588575.0000 18686 0.00%

### Conventions

Problem statistics: 352037 columns and 122283 rows.

231136 variables have integrality restrictions.

## Monday, February 9, 2009

### Set Covering Variant

*> I need to model a variant of the set covering problem:*

*> we need to cover at least p% of the points.*

The set covering problem looks like:

param a{I,J},binary;

var x{J},binary;

minimize cost: sum{j in J} x[j];

subject to covers{i in I}: sum{j in J} a[i,j]*x[j] >= 1;

We can formulate your problem with:

param a{I,J},binary;

param p >= 0, <= 1;

var x{J},binary;

var y{I} >= 0, <= 1;

minimize cost: sum{j in J} x[j];

subject to covers{i in I}: sum{j in J} a[i,j]*x[j] >= y[i];

subject to req: sum{i in I} y[i] >= p*card(I);

I believe y does not have to be binary and can be relaxed (between 0 and 1).

### Solver Status

- When building algorithms in GAMS. In this case the progress and flow of control is determined by the return codes from the SOLVE statement.
- When building user interfaces around a model. In this case we want to provide informative messages to the user.
- Or both (this was my case: solve is both embedded in an algorithm and I need to report useful messages to the user).

When relying on the Model Status and Solver Status in these cases, we assume that the solver is predictable. Here is a case where the GAMS/Cplex link is sloppy in reporting an iteration limit:

` S O L V E S U M M A R Y`

MODEL m OBJECTIVE z

TYPE MIP DIRECTION MINIMIZE

SOLVER CPLEX FROM LINE 27

**** SOLVER STATUS 4 TERMINATED BY SOLVER

**** MODEL STATUS 14 NO SOLUTION RETURNED

**** OBJECTIVE VALUE 0.0000

RESOURCE USAGE, LIMIT 226.084 1000.000

ITERATION COUNT, LIMIT 0 100000

ILOG CPLEX Dec 1, 2008 22.9.2 WIN 7311.8080 VIS x86/MS Windows

Cplex 11.2.0, GAMS Link 34

Cplex licensed for 1 use of lp, qp, mip and barrier, with 4 parallel threads.

MIP status(110): error termination, no integer solution

*** CPLEX Error 3019: Failure to solve MIP subproblem.

Error solving MIP subproblem.

Solution aborted due to iteration limit.

The correct solver status would be 2:

**** SOLVER STATUS 2 ITERATION INTERRUPT

Note also that the iteration count is not returned correctly. The link claims zero iterations have been performed. The correct number is 100000. This also means M.ITERUSD is wrong and can not be relied on when reporting back.

## Thursday, February 5, 2009

### floor in MINLP model

*> My MINLP model (in AMPL) does not solve:*

*> the relaxed NLP is declared infeasible.*

> The model contains the floor() function.

> The model contains the floor() function.

Unless you are using a global solver, MINLP solvers assume the relaxations (NLPs formed by fixing or relaxing the integer variables) are smooth. The floor function is discontinuous, so it is dangerous to use. There is a simple reformulation however, which moves the complexity from the NLP to the MIP part:

var x;

var flx integer;

flx ≤ x ≤ flx + 1

To get slightly more predictable behavior when x is (close to) integer, use:

flx ≤ x ≤ flx + 0.9999

In AMPL you may want to rewrite this as:

0 ≤ x − flx ≤ 0.9999

## Monday, February 2, 2009

### Social Golfer Problem in OML (2)

`Model[`

// N : number of golfers

// NG : number of groups

// T : number of rounds

// GS : group size (N/NG)

Parameters[Integers,N=16,NG=4,T=5,GS=4],

// Parameters[Integers,N=32,NG=8,T=10,GS=4],

Decisions[

Integers[0,NG-1],

Foreach[{i,N},{t,T},x[i,t]]

],

Constraints[

// form groups

Foreach[{g,NG},{t,T},Sum[{i,N},AsInt[x[i,t]==g]] == GS],

// golfer i and j meet at most one time

Foreach[{i,N}, {j,i+1,N}, Sum[{t,T},AsInt[x[i,t]==x[j,t]]] <= 1],

// fix first round

Foreach[{g,NG},{k,GS},x[GS*g+k,0]==g],

// fix first golfer

Foreach[{t,1,T},x[0,t]==0]

]

]