How to linearize

> I have this expression x₁⋅x₂⋅x₃, where all variables are binary. How can I linearize this?

The nonlinear equation z=x₁⋅x₂⋅x₃ with x_i ∈ {0,1} can be linearized as:

z ≤ x₁
z ≤ x₂
z ≤ x₃
z ≥ x₁+x₂+x₃−2

z can be continuous with z ∈ [0,1]. See also http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-binary-variables.html.

GAMS Solaris 64 bit problem

> 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.

Speeding Up Gams

The following calculations took more than 4 hours:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =
    sum(phzone$CountryLGSxTRCrops(faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr), 1);

scalar max;
max = smax((faoqcountry,faocrop,grid,irrg_key,lgsxtr), count(faoqcountry,faocrop,grid,irrg_key, lgsxtr));
display max;

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:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
scalar mx /0/;
count(faoqcountry,faocrop,grid,irrg_key,lgsxtr)=0;
loop((faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr)$CountryLGSxTRCrops(faoqcountry,faocrop,grid,phzone,irrg_key,lgsxtr),
    count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =  count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) + 1;
    mx = max(mx,count(faoqcountry,faocrop,grid,irrg_key,lgsxtr));
);

Explicit loops are expensive in GAMS but still much cheaper compared to dense processing. Even faster is:

parameter count(faoqcountry,faocrop,grid,irrg_key,lgsxtr);
parameter ReorderedCrops(faoqcountry,faocrop,grid,irrg_key,lgsxtr,phzone);
option ReorderedCrops < CountryLGSxTRCrops;

count(faoqcountry,faocrop,grid,irrg_key,lgsxtr) =
     sum(phzone$reorderedCrops(faoqcountry,faocrop,grid,irrg_key,lgsxtr,phzone), 1);

scalar max;
max = smax((faoqcountry,faocrop,grid,irrg_key,lgsxtr)$count(faoqcountry,faocrop,grid,irrg_key,lgsxtr),
             count(faoqcountry,faocrop,grid,irrg_key,lgsxtr));
display max;

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.

SQL2GMS tricks

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:

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

Note: we generate more meaningful set elements than 1, 2, 3,… by prefixing them with an ID string so they become faotqcntry_1, faotqcntry_2, faotqcntry_3,… This helps when debugging GAMS models.

See also [link].

Unjust benchmark: TSP MIP/Gurobi vs GA/Octave

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.

	MIP/Gurobi	GA/Octave
Source	eil51.gms	tsp_ga.m
Language	GAMS	Matlab
Solver	Gurobi	Octave
Options	threads 4 cuts 2 heuristics 0.2	defaults
Objective	426.00 (Proven optimal)	437.23
Time (seconds)	214.255	1939.2

Warning: to a large extent this is not a fair comparison.

TSP Tour from GAMS/Gurobi:

TSP Tour from TSP_GA:

Binary QP

I have been trying to solve an Binary Integer Program using AMPL/CPLEX
with the following objective function

minimize Total_Cost:
        sum{i in flights,k in gates} In[i,k] * Pnty * (Y[i,k]-X[i,k])^2 + sum
{i in flights, k in gates} ct*FTG[i,k]*(Y[i,k]-X[i,k])

where X[i,k] is the only decision variable.

I tried to solve it, but it gave me the following error:

CPLEX 11.2.0: 5 diagonal QP coefficients of the wrong
Variable Diagonal
_svar[8] -5
_svar[4] -2.5
_svar[3] -24.5
Diagonal QP Hessian has elements of the wrong sign.

When I tried to solve the problem using Minos (trial version) solver
for miniature problem, it solves it but the integrality of some
variables are ignored.

I am trying to linearize the above objective function and trying to
solve it through the CPLEX, but I am not successful yet. Can you
suggest for any solution or approach to solve it.

Assuming Y is a parameter and X is a binary variable (and all other identifiers are parameters), you can replace

(Y[i,k]-X[i,k])^2 = Y[i,k]-2*Y[i,k]*X[i,k]+X[i,k]^2 = Y[i,k]-2*Y[i,k]*X[i,k]+X[i,k]

which is linear in X. Note that x²=x when x is a binary variable.

24 cores

This [link] is to practice your French. Microsoft Solver Foundation on a machine with 24 cores to schedule the TechDays 2009 (Paris) conference.

GAMS/Gurobi link

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:

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

     MODEL   william             OBJECTIVE  cost
     TYPE    MIP                 DIRECTION  MINIMIZE
     SOLVER  GUROBI              FROM LINE  81

**** SOLVER STATUS     1 NORMAL COMPLETION        
**** MODEL STATUS      1 OPTIMAL                  
**** OBJECTIVE VALUE           988540.0000

---- EQU maxu  maximum generation level (1000mw)

                       LOWER          LEVEL          UPPER         MARGINAL

type-1.12pm-6am        -INF          -13.8000          .              .         
type-1.6am-9am         -INF           -8.0000          .              .         
type-1.9am-3pm         -INF          -13.0000          .              .         
type-1.3pm-6pm         -INF           -2.7500          .              .         
type-1.6pm-12pm        -INF          -12.7500          .              .         
type-2.12pm-6am        -INF           -0.4500          .              .         
type-2.6am-9am         -INF             .              .         -2100.0000  NOPT
type-2.9am-3pm         -INF             .              .         -4200.0000  NOPT
type-2.3pm-6pm         -INF             .              .         -2100.0000  NOPT
type-2.6pm-12pm        -INF             .              .         -4200.0000  NOPT
type-3.12pm-6am        -INF             .              .              .         
type-3.6am-9am         -INF             .              .              .         
type-3.9am-3pm         -INF             .              .              .         
type-3.3pm-6pm         -INF           -5.0000          .              .         
type-3.6pm-12pm        -INF             .              .              .         

**** REPORT SUMMARY :        4     NONOPT ( NOPT)

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).

Update: this is fixed.

glpk/ampl to GAMS conversion

> I don’t have glpk2gams, how do I convert an ampl or glpk model to gams?

1. Generate an MPS file from AMPL or GLPSOL:
1. Use the command line flag –wmps in glpsol
2. In AMPL you can use

ampl: option auxfiles rc;
ampl: write mxxx;

Here xxx can be replaced by a better name

 

2. Use mps2gms to translate this MPS file into a GAMS file. Note: mps2gms is part of the GAMS distribution so you already have it.

GAMS: set up a one-to-one mapping

In a model I was working on yesterday, I needed something like:

set i 'jobs' /start1*start10, job1*job100, end1*end10/;
set r 'resources' /r1*r10/;
set startmap(i,r) /
   start1.r1
   start2.r2
   start3.r3
   start4.r4
   start5.r5
   start6.r6
   start7.r7
   start8.r8
   start9.r9
   start10.r10
/;
set endmap(i,r) /
   end1.r1
   end2.r2
   end3.r3
   end4.r4
   end5.r5
   end6.r6
   end7.r7
   end8.r8
   end9.r9
   end10.r10
/;

This “diagonal mapping” can be written more compactly as using some new GAMS syntax:

set i 'jobs' /start1*start10, job1*job100, end1*end10/;
set r 'resources' /r1*r10/;
set startmap(i,r) /
   start1*start10:r1*r10
/;
set endmap(i,r) /
   end1*end10:r1*r10
/;