All-different

> Hi I am trying to solve a certain linear integer problem with lp_solve.
> I know that the solution set (24 variables) is (1,2...,24) and I am trying to figure out a way to input to lp_solve this array in order to get a feasible solution.

> In other words is it possible to put have a (or a set) of constraints saying that x1 ,x2,....,x24 belong to (1 ,2,...24) set but they are distinct (for instance x1=24, x2=3,x3=21 etc...).

This is called the ALL-DIFFERENT constraint. See here for an example.

One way of modeling this as a MIP is:

binary variable p(i,j) (i,j=1..24)

sum(i, p(i,j)) = 1;  for all j

sum(j, p(i,j)) = 1;  for all i

a(i) = i     (constant)

x(i) = sum(j, p(i,j)*a(j))

A tighter formulation would be:

sum(i, x(i)) = 0.5*n*(n+1), where n=24

sum(j in J, x(j)) ≥ 0.5*card(J)*(card(J)+1) where J is an element of the power set of {1,..,24} with card(J) < 24.

See this presentation.

Regular Expressions in GAMS/IDE

I had to add a ';' to each line in a file. In vi (unix editor) one would do:

1,$s/$/;/

In the GAMS IDE this works slightly different:

replace: .*$
by: $0;

As an aside: the equivalent of 1,$s/^/;/ is not obvious to me (there may be an issue here if we use ;$0 as replacement string).

AMPL NLP question

> I want to formulate some equations like ln(z(k)) <= D, 
> z(k) is a variable. Is there some function or keyword used as logarithm? 
> Or is there some method to calculate logarithm? Thanks 

You can use the log() function. Note that if D is a parameter the constraint can better be formulated as z[k] <= exp(D). This makes it a simple linear constraint that can be converted into a bound.

The difference can be illustrated with a small example:

ampl: var x := 1;
ampl: param D := 3;
ampl: maximize obj: x;
ampl: e: log(x) <= D;
ampl: solve;
MINOS 5.5: optimal solution found.
11 iterations, objective 20.08553692
Nonlin evals: constrs = 42, Jac = 41.
-----------------------------------------
ampl: var x := 1;
ampl: param D := 3;
ampl: maximize obj: x;
ampl: e: x <= exp(D);
ampl: solve;
MINOS 5.5: optimal solution found.
1 iterations, objective 20.08553692

AMPL has a very good presolver, but not so good this reformulation is applied automatically.

Gurobi vs Cplex benchmarks

> What are your first impressions of Gurobi? How does it compare to CPLEX for MIP? Have you done any benchmark analysis?

Here are some benchmarks: http://plato.asu.edu/ftp/milpc.html

However I would always urge you to try out solvers on your particular models as they may behave very differently than the benchmark models. Drop me a note if you need me to run some models or need access to evaluation systems.

Cplex log: number of integer variables

Sometimes progress logs pose some interesting questions. In the Cplex log below we solve a problem with 2016 binary variables (after presolving). However the progress log seems to indicate that when we start 2036 of them are integer infeasible.

Reading data...
Starting Cplex...
Tried aggregator 1 time.
MIP Presolve eliminated 45077 rows and 15073 columns.
MIP Presolve modified 504 coefficients.
Aggregator did 1008 substitutions.
Reduced MIP has 73884 rows, 26656 columns, and 201936 nonzeros.
Reduced MIP has 2016 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time =    0.45 sec.
Clique table members: 900.
MIP emphasis: hidden feasible solutions.
MIP search method: dynamic search.
Parallel mode: opportunistic, using up to 4 threads.
Parallel mode: opportunistic, using up to 4 threads for concurrent optimization.
Using devex.
Using devex.
 
Total real time on 4 threads =   16.18 sec.
Root relaxation solution time =   16.19 sec.
 
        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer     Best Node    ItCnt     Gap
 
      0     0        0.0000  2036                      0.0000       20         
 

Adding redundant constraint to MIP

Whenever a MIP model is very slow I try to invent "redundant" constraints that can be added to the model. Here is an example with Cplex. We report number of nodes needed to prove optimality on some small examples of a (proprietary) scheduling model:

case	with	without
8-3-2	23	34
8-6-2	0	7290
9-3-3	0	81
9-4-3	0	180
9-5-3	74	87738

As can be seen the version with the redundant equation added has much better performance. Furthermore it is noted that the duality gap was much smaller.

Scheduling: limits on number of periods active.

> How to model that an employe must work at least as 5 successive periods  and at most as 8 successive periods with mathprog.

There are a few ways to model this, and without knowing more about the model it is difficult to give definite advise, but here is a suggestion of a possible formulation taken from a machine scheduling model I developed a little while ago:

param MAXT := 25;
param MIN := 5;
param MAX := 8;
set T := {1..MAXT};
 
var x{T}, binary;
var SwitchOn{T}, binary;
var SwitchOff{T}, binary;
 
 
 switch1{t in T: t > 1}: x[t]-x[t-1] = SwitchOn[t]-SwitchOff[t];
 switch2{t in T}: SwitchOn[t]+SwitchOff[t] <= 1;
 
 init1: x[1] = 0;
 init2: SwitchOn[1] = 0;
 init3: SwitchOff[1] = 0;
 
 
 minimum{t in T}: sum{i in 1..MIN:t+i-1<=MAXT} x[t+i-1] >= MIN*SwitchOn[t];
 maximum{t in T:t+MAX<=MAXT}: x[t+MAX] <= 1-SwitchOn[t];

I ignored the details about what to do with the boundaries. The last condition adds also an implied restriction on when the machine can be turned on again (after 8 periods), so it may be better to model this as:

maximum{t in T:t+MAX<=MAXT}: sum{i in 1..MAX+1:t+i-1<=MAXT} x[t+i-1] <= MAX + 1 - SwitchOn[t];

For multiple persons j we would need x[j,t], SwitchOn[j,t], SwitchOff[j,t]. A different way would be to organize x[j,t,len] where t is the start time and len is the length of occupation (len=5,6,7,8). This approach has more integer variables but may solve fast depending on the situation.

Excel Solver Sensitivity Analysis

> How do I programmatically tell Solver to generate a sensitivity report

> after running SolverSolve ?

This can be done through the SolverFinish function.

Updated lp2gams

Updated lp2gams to distinguish names with different casing. In lpsolve a variable x is different from a variable X while in GAMS they are the same. We added some code to make sure the generated GAMS names are unique. If needed we add an underscore to the name.

Cplex barrier results

With some simulated data for the problem described in this post, we achieved some very good results with the Cplex barrier method. Below are timings in seconds with different Cplex algorithms

	rows/cols/nz	dual	primal	barrier
model 1	1e5/2e5/1.4e6	1247	857	8.7
model 2	1.1e5/2.1e5/6.1e5	195	1223	7.3

The GAMS models look like:

41  variable z,a1(j),a2(j),s1(i),s2(i),y(k);
42  positive variable s1,s2;
43  equation obj,e1(i),e2(k),e3(i);
44
45  option lp=cplex;
46  option iterlim=1000000;
47  option reslim=10000;
48
49  obj.. z =e= sum(i, s1(i)+s2(i));
50  e1(i).. sum(j, xx(i,j)*a1(j)) + s1(i) - s2(i) =e= b(i);
51  model m1 /obj,e1/;
52  solve m1 minimizing z using lp;
53
54  e2(k).. y(k) =e= sum(j, x(k,j)*a2(j));
55  e3(i).. sum(map(i,k), y(k)) + s1(i) - s2(i) =e= b(i);
56  model m2 /obj,e2,e3/;
57  *solve m2 minimizing z using lp;

Note that the SUM in line 55 is actually not a real summation. It just selects a k belonging to a particular i. A different, maybe better, way to write this would be:

equation e3(i,k);

e3(map(i,k)).. y(k) + s1(i) - s2(i) =e= b(i);