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

Question

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]
g = o.add_var(name='max_i', var_type=CONTINUOUS)
h = o.add_var(name='min_i', var_type=CONTINUOUS)

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

Answer 1

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;

Answer 2

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