# 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, 2021 at 23:34
• Not sure, but nz count reduction was the motivation. Aug 6, 2021 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, 2021 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, 2021 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, 2021 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, 2021 at 21:20
• It shows modell has 282 vars, 154 constraints and 13960 nzs, Aug 5, 2021 at 21:26
• That is very small. (With #nz's a bit too big for my taste) Aug 5, 2021 at 21:27