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, April 30, 2009

I´m making my final project working with AMPL and nobody in my university can help me. […] The question is the following: I must take a decision over a variable, DECISION[n]. This variable can be a number between 0 to 64. I want to obtain another variable, RESULT[n], that if the first variable is Zero, the second one is zero too. In the other cases, must be 1. In C language it could be: if DECISION[n] =0 then RESULT[n]=0 (else RESULT[n]=1). i´m working with CPLEX.

E.g. using Cplex indicator constraints:

EQ{n in N}: RESULT[n]=0 ==> DECISION[n]=0 else DECISION[n]>=0.0001;

You really should get advise from your professor. He/she gets paid to help you with this.

Wednesday, April 29, 2009

I was involved in a small project where we converted a spreadsheet model using Solver to GAMS/CONOPT. I was unable to reproduce the spreadsheet results, and it turns out the spreadsheet had some errors in the formula’s. I have converted quite a number of Excel models, and almost all had some problems that were uncovered by moving the model to a modeling language environment. See table 1 in http://panko.shidler.hawaii.edu/SSR/Mypapers/whatknow.htm for more evidence of how difficult it is to create error-free spreadsheet models.

PS: after fixing the spreadsheet errors we were able to get similar results. This is actually not a bad approach to verify and repair spreadsheet models.

Hello, I have to deal with an objective function of the following form: min [C1 * ceiling (X/m) plus C2 * (X-ceiling(X/m))]; where C1, C2, m are parameters and X is the decision variable. Ceiling is a mathematical operator that rounds a number up, to the nearest integer (see Excel's definition). My first question is: is the Ceiling operator defined in OML? My second question: is there a solver within the Microsoft Solver Foundation that can handle this type of objective function? Thank you.

I believe you can do the following in OML:

Let Y be an integer variable with the constraint

X/m ≤ Y ≤ 1+X/m

then minimize

C1 * Y + C2 * (X-Y)

Hopefully the rest of the model is linear. In that case this is a much better formulation as it gets rid of a non-smooth function: solvers often have troubles with that.

Actually I would have expected problems in the big equality condition and not in the Unequal constraint.

Update: I understand now the behavior a little bit more. Without the Unequal the MIP solver Gurobi is called. This produces the correct conclusion: Feasible. With the Unequal statement enabled the CSP solver is used and that solver has problems with this model. Actually Gurobi uses floating point instead of rational arithmetic. So it is somewhat surprising that Gurobi is more accurate on this model than the CSP solver.

Update2: see expert explanation in comment. CSP wants all intermediate results (I think all terms used in any expression) also to fit in an integer.

The numbers in the addition equation are really become too big for comfort. It is surprising that some solvers such as Gurobi can deliver a good solution for this model, but as expected some others can not. The poor scaling can create severe havoc. The model can be downloaded from here.

Saturday, April 25, 2009

A different formulation can be developed that circumvents this problem. We can use the addition algorithm taught to us in primary school: first add the right most digits, and save a possible carry. Then add up the digits in the second rightmost column, including the previous carry, etc. This scheme can actually be implemented in constraints as follows:

This solution is correct. At this moment I have no explanation why the GAMS version is doing worse. There may be issues with compilers, GAMS handling of the objective or other possible reasons. There is a difference in the log indicating the scaling is done differently. I believe the standalone version is 4.34 (Dec 2008) compared to 4.32 (built Feb 12 2009) for the GAMS version, so the GAMS system was little bit outdated. Also note that I used a 64 bit GAMS system and a 32 bit GLPSOL system. It remains very treacherous to just compare results from a numerical algorithm, but in this case we expect things to be more equal than they are now, and it may be possible we stumbled upon a bug.

with the variables yi∈{0,1,...,9}. This is not sufficient as a solution yi=0 would be allowed. In addition we need to add an all-different constraint, and we may also need to forbid leading zero’s.

Forbidding leading zero’s is easy. We add the bounds: ys≥1 and ym≥1. The other constraint is somewhat more complicated. The uniqueness constraint can be implemented using an extra binary variable xi,j with form a permutation matrix.

We need all 9 y's so we introduce ydummy1 and ydummy2. This now results in the MIP model:

Note that GAMS does not report zero’s. So you will need to deduct that the missing y(‘O’) is zero. This model solves quickly with all GAMS solvers except BDMLP (the LP seems to cycle on one of the nodes).

The objective function is somewhat interesting: the objective of the relaxed problem is always identical to the objective of the MIP. We can see this from the following:

As a result the MIP solver can stop as soon as it finds the first integer solution: it is immediately optimal. Of course a simpler way of achieving the same would be to use the following objective:

The graph in this model is somewhat misleading as the coordinates are in DD.MM format where DD are degrees and MM are minutes. A more realistic graph is created by using decimal degrees: DD + MM/60. Note: the optimal tour is calculated correctly as the distances were ok. Just the display of cities was slightly inaccurate. See below:

Thursday, April 16, 2009

> When I run GAMS on a Solaris 64 bit SPARC system I encounter the following problem: > > ld.so.1: gmsco3ux.out: fatal: libfui.so.1: open failed: No such file or directory

This means you are missing Fortran run time support. A client of mine had the same problem when I sold him a GAMS system for Solaris 64-bit SPARC. In this case you try to run the CONOPT solver, which is written in Fortran. It is dynamically linked and therefore needs some run time libraries. If you don't have a Fortran compiler installed on your system you may get this error. Software suppliers are allowed to provide these files. The GAMS web site actually does provide a download for the 32-bit runtime here. Unfortunately for 64 bit no such download is provided. You can contact GAMS support to request the 64 bit Fortran run time libraries.

The parameter CountryLGSxTRCrops is very large (but sparse): 2.5 million elements. The reason for this fragment to be so slow is two-fold:

In the assignment to parameter count we are doing things out of order. GAMS prefers to loop over sets at the end of the index list, but here we loop over an index in the middle (phzone). Sometimes GAMS is able to reorder things automatically to do things faster, but in this case apparently not.

smax is executed dense, as zeroes can be significant. In GAMS zero and “does not exist” is the same. This can be exploited when sparse processing is used. Here GAMS is conservative and reverts to dense processing as the implicit zeros may be important in calculating the correct SMAX value.

This a direct rewrite into an explicit loop that executes in 30 seconds:

Here we use the (undocumented) < option to reorder the big parameter such that the phzone index is last. This is the representation that GAMS likes most when looping over phzone. Now the calculation of count is very fast. Next we force the calculation of max to be performed over the nonzero elements in parameter count only. This formulation takes less than 5 seconds. We use some more memory by duplicating CountryLGSxTRCrops.

Sunday, April 12, 2009

Using SQL2GMS to read a large set of CSV files has the big advantage to be able to repair things on the fly. We had a few NULL values in one of the CSV files and I fixed this by skipping these records:

Q39=select 'faotqcntry_'&intFAOLVSTQ_CNTRY,'faolivestp_'&intFAOLVST2_CODE,[dblFAOLVSTPRICE_$S/MT] from LivestockPrices.csv where not isnull([dblFAOLVSTPRICE_$S/MT])

Then I got the request to use a value of 99999 for those records. This can be easily done by:

Wednesday, April 8, 2009

I was asked to have a look at the Genetic Algorithm from this site: http://www.mathworks.com/matlabcentral/fileexchange/13680. Just to get an idea I benchmark here the TSP model eil51 from TSPLIB against GAMS/Gurobi. I used Octave instead of Matlab. The solution TSP tour is quite different.

Sunday, April 5, 2009

During a demo of GAMS/Gurobi we found a minor error. It shows NOPT's for equations in optimal MIP models. No show-stopper: you should just ignore these messages for now. This is with model magic.gms from the model library:

Note: these NOPTS are used to mark rows and columns with the wrong sign for the marginals in LP's. In this case the signs for these row marginals (i.e. duals) look just fine. In GAMS, after a MIP is solved, an LP is formed by fixing the integer variables. The marginals reported in the listing file is for this fixed LP problem. In general, NOPT flags should never appear in a model that is declared optimal: they are reserved for models that are intermediate non-optimal, a status that can happen e.g. when hitting a time or iteration limit before the model was optimal. If we analyze this a little bit further, we note that for an optimal minimization problem all marginals corresponding to non-basic rows (and columns) at upper bound should be negative or EPS. This is the case for these rows, and the NOPT flag is therefore inaccurate. Note that a row or column can be basic while at bound (this is also known as degeneracy). We see a few of those degenerate rows also (they are reported correctly).