Here is a toy model that illustrates that certain models are poorly suited for solving by MIP solvers.
Sorting data is a rather elementary operation. Surprisingly, GAMS has nothing built in for this. There is an external tool for sorting (gdxrank), but in some environments (e.g. NEOS) this is not allowed. There are some slower GAMS coding tricks that can help. Of course we can also write a simple MIP model that performs sorting:
| $ontext |
This is quite slow however. n=100 is no problem, but n>200 will require significant time (>1000 seconds). Constraints assign1 and assign2 form an assignment problem, which is easy to solve. The simple extension will cause relaxation based MIP solvers to slow down enormously. Obvious question: are there better formulations? I would guess constraint programming solvers would do better on this.
For some instances Gurobi is doing very good, probably because of its excellent heuristics. But this is not the case for all instances: in some cases it is just as slow as other solvers.
The above model contains two noteworthy modeling tricks:
- a dummy objective so that we stop after the first integer solution is found
- the equation sort is indexed over i+1 so that the last one is not generated: this equation has card(i)−1 members.
I am glad that I never tried to learn GAMS :) can't understand anything from the code.
ReplyDeleteThe second issue is probably a result of the first one.
ReplyDeleteErwin,
ReplyDeleteI don't know much about GAMS, but I think I understand your model.
I'd like to propose a "Sorting by LP" formulation for the same problem. It uses the objective function in a simple but effective way to obtain the correct assignment.
Since the model basically forms a network model (leaving out the helper variables y which only facilitate the output), all continous variables take integer values in the optimal solution.
Using this formulation, CPLEX 12.3 solves a problem with n = 500 within 5 seconds using the barrier solver using 2 threads.
This is the model in AMPL / GMPL:
param n := 500 integer;
set INDEX := 1..n;
param value{INDEX} := Uniform01();
var assign{INDEX,INDEX} >= 0, <=1;
var y {INDEX} ;
maximize obj:
sum{i in INDEX, j in INDEX} assign[i,j] * i * value[j];
subject to assign1{i in INDEX}:
sum {j in INDEX} assign[i,j] = 1;
subject to assign2{j in INDEX}:
sum {i in INDEX} assign[i,j] = 1;
subject to y_def{i in INDEX}:
y[i] = sum {j in INDEX} assign[i,j] * j;
solve;
display y;
Great! May be divide value[j] in obj by max of values, so it also works for numbers > 1.
ReplyDeleteActually, I was shortly thinking about the same issue, but the range of the values does not matter (even if the values are negative). It is only important that the factor, in this case j, is strictly ordered.
ReplyDeleteSorry, in the previous post I meant the factor i which occurs in the objective function. You
ReplyDeletecould replace the i with any factor factor[i] as long as factor[i] > factor[i-1] for all i and factor[i] >= 0.
The solver will always allocate the j with the maximum value[j] to the biggest factor since this has the greatest impact on the objective function (and so on...)
Got it!
ReplyDeleteWhat would happen if you replaced the sort() constraint by y(i) >= y(j) for all i > j? Of course it would make the model larger, but would it make the model easier to solve? Another idea: create an additional set of binary variables that are 1 if y(i+1) >= y(i) and 0 otherwise, then maximize the sum of these binary variables.
ReplyDeleteI have made some further tests with CPLEX 12.3, using a model with n = 1000, that is, we have 1 Million variables (+ the 1000 helper variables which are typically eliminated in the presolve phase), on an Intel dual core machine. I have tested with 3 different LP algorithms:
ReplyDeletePrimal Simplex (CPLEX's automatic choice): 92 sec.
Barrier: 22 sec.
Network Simplex: 4 sec.
I append the model again, having added some comments and the variable Sorted which has the same meaning as Erwins variable y.
#=============================DECLARATION
param n := 500 integer;
set INDEX := 1..n;
param values{INDEX} := Uniform01(); # contains the unsorted values
var Assign{INDEX,INDEX} >= 0; # assigns the (sorted) index i to the value with index j
var Sorted{INDEX}; # helper variable: will contain the sorted values.
#==============================MODEL
# this objective does the whole trick:
# by using i as a factor, the j with the greatest value will be assigned to the
# greatest index n, the second greatest value to index n-1 and so on
maximize obj: sum{i in INDEX, j in INDEX} Assign[i,j] * i * values[j];
# the assignment constraints
subject to assign1{i in INDEX}: sum {j in INDEX} Assign[i,j] = 1;
subject to assign2{j in INDEX}: sum {i in INDEX} Assign[i,j] = 1;
# this constraint only computes the sorted array
subject to sorted_def{i in INDEX}: Sorted[i] = sum {j in INDEX} Assign[i,j] * values[j];
#==============================SOLVE AND DISPLAY
solve;
display values, Sorted;
replace the sort() constraint by y(i) >= y(j) for all i > j.
ReplyDeleteIn GAMS:
sort(i,j)$(ord(i)>ord(j)).. y(i) =g= y(j);
I don't think this will cut off any fractional solutions, but may be it can help Gurobi creating better cuts.
Here are the results with Gurobi:
n=100, old sort: 6.8 sec, new sort: 37.3 sec
n=150, old sort: 22.0 sec, new sort: >1000 sec
So this does not seem to help (although 2 models is hardly proof).