How to solve MILP efficiently using MIP python? Any other solvers to quickly solve where I have only integer and binary variables alone?

from mip import *

q=[[4,5,0,2],[3,2,1,0],[0,1,4,5],[3,2,0,1],[2,1,3,2],[4,3,2,1],[3,4,2,0],[5,0,3,1],[0,4,3,3],[4,0,4,1]]
s=[11,1,1,1,1,1,1,1,1,1]
r=[6,8,4,7]
a=range(0,10)
b =range(0,20)
c = range(0,4)
d =[12,9,6,6]
o = Model()
x= [[o.add_var('x({},{})'.format(i, j), var_type=BINARY) for j in b] for i in a]
y = [[o.add_var('y({},{})'.format(k, j), var_type=BINARY) for j in b] for k in c]

for i in a:
o += xsum(x[i][j] for j in b) == s[i]
for j in b:
o += xsum(x[i][j] for i in a) == 1
for k in c:
o += xsum(y[k][j] for j in b) == d[k]
for j in b:
for k in c:
o += r[k] *(xsum(y[k][jj] for jj in b if jj <= j)) >= (xsum(( xsum (q[i][k]* x[i][jj] for i in a )for jj in b if jj <=j)))
for j in b:
o += (xsum( xsum(r[k] * y[k][jj] for k in c) for jj in b if jj <= j)) - (xsum(xsum( xsum (q[i][k]* x[i][jj] for i in a) for k in c) for jj in b if jj<=j)) <= g
for j in b:
o += (xsum( xsum(r[k] * y[k][jj] for k in c) for jj in b if jj <= j)) - (xsum(xsum( xsum (q[i][k]* x[i][jj] for i in a) for k in c) for jj in b if jj<=j)) >= h
o.objective = minimize(g-h)
o.optimize()
print("Solution with {} found.".format(o.objective_value))
for x in o.vars:
print('{} : {}'.format(x.name, x.x))


"""I have tried to solve the above MILP using python mip and tried cbc, gurobi solvers, i can able to solve it for small problem , lets say when s=[20,0,0,0,0,0,0,0,0,0] and d = [14,13,0,6], and also when length of s is less than 5, but for the above inputs, it taking too much time (more than 8 hours and i can't wait).

In the above test case, I have 282 variables and 154 constraints only and I have seen that number of variables and number of constraints alone doesn't contribute to how much time, the solver takes to solve. I have 10k variables and 7k constraints in my main problem to solve for. I have tried constrained programming by limiting the values of variables by lb and ub, it didnt help. How can I solve it better? Help would be appreciated. ( g and h will be integers in all case, so I tried it by initializing as integer variables and giving close lb and ub, It also didn't help)"""

Pushing @Erwin's idea further, I got faster solve times by instead introducing variables ry and qx:

for j in b:
for k in c:
o += ry[j][k] == r[k] *(xsum(y[k][jj] for jj in b if jj <= j))
for j in b:
for k in c:
o += qx[j][k] == xsum(( xsum (q[i][k]* x[i][jj] for i in a )for jj in b if jj <=j))
for j in b:
for k in c:
o += ry[j][k] >= qx[j][k]
for j in b:
o += (xsum(ry[j][k] for k in c)) - (xsum(qx[j][k] for k in c)) <= g
for j in b:
o += (xsum(ry[j][k] for k in c)) - (xsum(qx[j][k] for k in c)) >= h


The optimal objective value turns out to be $$3$$, with $$(g,h)=(15,12)$$.

As requested, here is the full SAS code, with variables ry[k,j] and qx[k,j]:

proc optmodel;
set ASET = 0..9;
set BSET = 0..19;
set CSET = 0..3;
num q {ASET, CSET} = [
4,5,0,2,
3,2,1,0,
0,1,4,5,
3,2,0,1,
2,1,3,2,
4,3,2,1,
3,4,2,0,
5,0,3,1,
0,4,3,3,
4,0,4,1
];
num s {ASET} = [11,1,1,1,1,1,1,1,1,1];
num r {CSET} = [6,8,4,7];
num d {CSET} = [12,9,6,6];

var x {ASET, BSET} binary;
var y {CSET, BSET} binary;
var g;
var h;

con Con1 {i in ASET}:
sum {j in BSET} x[i,j] = s[i];
con Con2 {j in BSET}:
sum {i in ASET} x[i,j] = 1;
con Con3 {k in CSET}:
sum {j in BSET} y[k,j] = d[k];

var ry {CSET, BSET};
con ryCon {k in CSET, j in BSET}:
ry[k,j] = r[k] * sum {jj in BSET: jj <= j} y[k,jj];
var qx {CSET, BSET};
con qxCon {k in CSET, j in BSET}:
qx[k,j] = sum {i in ASET, jj in BSET: jj <= j} q[i,k] * x[i,jj];
con Con4 {k in CSET, j in BSET}:
ry[k,j] >= qx[k,j];
con Con5 {j in BSET}:
sum {k in CSET} (ry[k,j] - qx[k,j]) <= g;
con Con6 {j in BSET}:
sum {k in CSET} (ry[k,j] - qx[k,j]) >= h;

minimize Range = g - h;
solve;

print x;
print y;
print g h;
quit;


And here is the resulting output:

x
0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
0  0  1  0  1  1  1  1  0  1  1  1  1  0  1  1  0  0  0  0  0
1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1
2  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0
4  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0
5  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0
6  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0
7  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
8  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0
9  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

y
0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
0  1  1  0  1  0  1  0  1  1  0  1  0  0  1  1  0  1  0  1  1
1  1  0  0  1  1  1  0  0  0  1  1  0  1  0  1  1  0  0  0  0
2  1  1  0  0  0  0  1  0  0  0  0  1  0  1  0  0  1  0  0  0
3  1  0  1  0  0  0  1  0  1  0  0  1  0  0  0  0  0  1  0  0

g  h
15 12


An alternative reformulation for your last three constraints is as follows:

   var ry_minus_qx {CSET, BSET} >= 0;
con Con4 {k in CSET, j in BSET}:
ry_minus_qx[k,j] = r[k] * y[k,j] - sum {i in ASET} q[i,k] * x[i,j] + (if j-1 in BSET then ry_minus_qx[k,j-1]);
con Con5 {j in BSET}:
sum {k in CSET} ry_minus_qx[k,j] <= g;
con Con6 {j in BSET}:
sum {k in CSET} ry_minus_qx[k,j] >= h;

• Just reduction in nz count? Or are there also other things going on. Aug 5 at 23:34
• Not sure, but nz count reduction was the motivation. Aug 6 at 0:16
• Hi, @RobPratt, thanks for the resopnse. I have tried the above code snippet, it didnt work for me and i think ry, qx index are interchanged and in the final two constraints, it is ry instead of r. SO, i modified to following and it also didn't work for me Aug 6 at 4:06
• ry = [[o.add_var('y({},{})'.format(k, j), var_type=CONTINUOUS) for j in b] for k in c] qx = [[o.add_var('y({},{})'.format(k, j), var_type=CONTINUOUS) for j in b] for k in c] for j in b: for k in c: o += ry[k][j] >= qx[k][j] for j in b: o += (xsum(ry[k][j] for k in c)) - (xsum(qx[k][j] for k in c)) <= g for j in b: o += (xsum(ry[k][j] for k in c)) - (xsum(qx[k][j] for k in c)) >= h  Aug 6 at 4:07
• Can you help me with how you solved, how you declared qx, ry and which solver you used, whats the solver time? Aug 6 at 4:10

How big is the model? Is the time spent in the solver or in generating the model? Can you should a bit of the solver log file (say the beginning and the end)?

I see one thing I don't like:

for j in b:
o += (xsum( xsum(r[k] * y[k][jj] for k in c) for jj in b if jj <= j)) - (xsum(xsum( xsum (q[i][k]* x[i][jj] for i in a) for k in c) for jj in b if jj<=j)) <= g
for j in b:
o += (xsum( xsum(r[k] * y[k][jj] for k in c) for jj in b if jj <= j)) - (xsum(xsum( xsum (q[i][k]* x[i][jj] for i in a) for k in c) for jj in b if jj<=j)) >= h


This essentially duplicates a sizable expression. I would use:

for j in b:
o += z[j] == (xsum( xsum(r[k] * y[k][jj] for k in c) for jj in b if jj <= j)) - (xsum(xsum( xsum (q[i][k]* x[i][jj] for i in a) for k in c) for jj in b if jj<=j))
o += z[j] <= g
o += z[j] >= h


(More equations and variables but fewer nonzero elements).

• Thanks Erwin for responding. I have tried already with sizable expression for constraints. The above code is my model for the test case, the time is being spent on the solver. Aug 5 at 21:20
• It shows modell has 282 vars, 154 constraints and 13960 nzs, Aug 5 at 21:26
• That is very small. (With #nz's a bit too big for my taste) Aug 5 at 21:27