Expected ILP solving time and how to improve speed

Question

I am trying to solve this ILP with all binary variables:

$$\min \sum_{1\leq i<j\leq n-1-h} t_{i,j}$$ $$\sum_{i=1}^{n-1-h} a_{k,i} = \lfloor (n-1)/2\rfloor \qquad \text{for }k\in[h];$$ $$a_{k,i} + a_{k,j} \leq 2 d_{k,i,j}\qquad \text{for }1\leq i<j\leq n-1-h,\ k\in[h];$$ $$\sum_{k=1}^h d_{k,i,j} \leq h - 1 + t_{i,j}\quad \text{for }1\leq i<j\leq n-1-h$$

thanks to an answer to my problem here, for $n=53$ and $h=13$. If I compute the minimum in the original problem by means of many random cases I have a minimum value of $113$, therefore the actual minimum will be less or equal than that.

I am using the free program LPSolve IDE 5.5.2.11 and the full problem for my $n=53$ and $h=13$ case is here (too long to copy within the question). I have encoded multiple variable indexes in one unique integer index for the software model.

There are 10387 constraints and 10881 variables.

The partial result after a couple of minutes was this:

Model name:  'LPSolver' - run #1
Objective:   Minimize(R0)

SUBMITTED
Model size:    10387 constraints,   10881 variables,        39780 non-zeros.
Sets:                                   0 GUB,                  0 SOS.

Using DUAL simplex for phase 1 and PRIMAL simplex for phase 2.
The primal and dual simplex pricing strategy set to 'Devex'.

Relaxed solution                   0 after      10774 iter is B&B base.

Feasible solution                320 after      22515 iter,      7838 nodes (gap 32000.0%)
Improved solution                319 after      23350 iter,      8716 nodes (gap 31900.0%)
Improved solution                318 after      24221 iter,      9633 nodes (gap 31800.0%)
Improved solution                317 after      25192 iter,     10569 nodes (gap 31700.0%)
Improved solution                316 after      26182 iter,     11556 nodes (gap 31600.0%)
Improved solution                315 after      27604 iter,     12611 nodes (gap 31500.0%)
Improved solution                314 after      28657 iter,     13614 nodes (gap 31400.0%)
Improved solution                313 after      29865 iter,     14645 nodes (gap 31300.0%)

Now the software has done about 4,000,000 iterations in about two hours without further progress.

I have very little knowledge of linear programming, so I am asking:

Is it possible to make that ILP problem faster by modifying it, use another software or change solver parameters?

Can we say anything about the expected completion time?

Answer 1

Is it possible to make that ILP problem faster by modifying it, use another software or change solver parameters?

Using other software: Definitely, as Rob Pratt's answer indicates. (I got an incumbent solution with value 67 in about 30 seconds with CPLEX, after which progress stalled.)

Changing solver parameters: This depends to some extent on the solver (different solvers have different parameters). For example, I set the "MIP emphasis" parameter in CPLEX to 5 (heavy emphasis on heuristics) to encourage it to get a "good" solution quickly.

Modifying the model: This is the tricky part. One problem with your model is that it has a "weak" formulation, meaning that relaxing integrality and solving it as an LP produces an answer with a very low objective value (0).* I tried switching the MIP emphasis in CPLEX to improvement of the lower bound, but after five minutes (the limit of my patience) the lower bound was still 0. So a stronger (tighter) formulation would definitely help ... if you could find one.

An obvious weakness is the constraint $a_{k,i} + a_{k,j} \leq 2 d_{k,i,j},$ which is presumably intended to say that $d_{k,i,j}$ must be 1 if either $a_{k,i}$ or $a_{k,j}$ is 1. You can get the same restriction using two constraints ($a_{k,i}\le d_{k,i,j}$ and $a_{k,j} \le d_{k,i,j}$ at the cost of a modest increase in model size. Using separate constraints is somewhat stronger. For instance, if in an LP relaxation the solver sets $a_{k,i}=1$ and $a_{k,j}=0,$ the revised formulation forces $d_{k,i,j}=1$ whereas the original formulation only forces $d_{k,i,j}=0.5.$ Unfortunately, in my tests, that change did not accomplish anything useful by itself.

Another issue with the model is symmetry. Branch and bound solvers can bog down when equivalent solutions exist in multiple branches of the search tree. In your case, permuting the rows of the $A$ matrix in any solution produces another solution with identical objective value.

Can we say anything about the expected completion time?

You will get a provable optimum within your lifetime ... if you are immortal (and assuming your computer does not run out of power or memory). Seriously, predicting solution times for integer programs is somewhere between impossible and not even that easy.

* Unless the actual optimal value is 0, in which case the problem is not a weak formulation but rather difficulty improving the primal bound (best feasible solution).

Answer 2

Here are the $a_{k,i}$ values for a feasible solution with objective value $61$, obtained in a few minutes with a different solver. Please verify.

1 0 1 1 1 1 0 1 1 1 0 1 1 0 1 0 1 0 1 1 0 1 1 1 1 1 1 0 0 0 1 1 0 0 1 1 0 1 1 
1 1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 0 0 1 0 0 0 
1 1 1 0 0 1 0 1 1 1 1 0 0 1 0 1 1 1 0 0 0 1 0 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 0 
0 0 1 1 1 0 1 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 1 0 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 
1 1 1 1 1 0 1 0 0 0 1 0 0 1 1 1 1 0 1 1 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 0 1 
1 1 0 1 1 1 1 1 0 1 0 1 0 1 1 0 1 1 1 0 0 0 1 0 1 1 1 0 1 1 0 1 1 1 1 0 1 0 1 
0 1 0 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 0 1 1 1 1 1 1 0 0 1 1 1 1 1 0 1 0 0 0 1 1 
0 1 1 0 1 1 0 0 1 0 1 1 1 1 0 1 0 1 1 1 1 1 1 1 0 0 1 0 0 1 0 1 1 1 0 1 1 1 1 
1 1 1 0 0 1 1 1 1 0 1 1 1 1 1 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 1 1 1 1 0 1 1 1 1 
0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 0 1 1 1 0 0 0 0 1 1 1 0 1 1 0 
1 1 0 1 0 0 0 0 1 1 1 0 1 0 1 1 0 1 1 1 1 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 
0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 0 1 1 0 1 1 0 1 0 1 1 1 1 1 1 1 0 0 1 0 1 0 1 
1 0 1 1 1 1 1 0 1 0 1 1 1 0 1 1 1 1 1 0 1 0 1 1 1 0 0 1 1 1 1 0 0 0 1 1 0 1 0 

---

Here's a better solution with objective value $57$, obtained by relaxing the first constraint from $=$ to $\ge$. Each row still has exactly $26$ ones.

1 1 0 0 1 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 0 1 1 0 1 1 0 1 0 1 0 1 1 1 0 0 1 0 1
1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 1 0 1 1 1 1 0 1 1 0 0 0 1 1 0 1 0 1 1 1 0 0 1 1
0 1 1 0 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 0 0 1 1 0 1 1 0 0
1 0 1 1 1 1 0 0 1 0 1 1 1 0 0 1 1 1 0 1 1 1 1 1 1 1 1 1 0 0 1 1 0 1 1 0 0 1 0
1 0 0 1 1 1 0 1 0 1 1 1 1 1 1 1 0 1 1 0 0 1 1 0 1 1 1 0 1 1 0 1 1 0 1 1 0 0 1
1 1 1 1 0 1 1 1 0 1 1 0 1 1 1 0 1 1 0 0 1 1 0 0 1 1 1 0 1 0 1 0 1 0 1 0 1 1 1
1 1 1 1 0 0 1 1 1 1 0 0 1 1 0 1 0 1 1 1 1 1 0 1 0 1 1 1 0 1 0 1 1 0 0 1 1 1 0
0 0 1 1 1 0 1 1 1 0 1 1 0 1 1 1 1 0 0 0 1 0 1 1 0 1 1 0 0 1 1 0 1 1 1 1 1 1 1
0 1 1 1 1 1 0 1 1 0 0 1 0 0 1 1 0 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 1 1 1 0 1
1 1 0 1 1 0 1 0 0 1 1 1 1 1 1 0 1 0 1 1 0 0 0 0 0 0 1 1 1 1 1 1 0 1 1 1 1 1 1
1 0 1 1 1 1 1 0 0 1 0 1 1 1 0 1 0 1 0 0 1 1 1 1 1 0 0 0 1 1 1 1 0 1 1 1 1 1 0
1 0 1 0 0 1 1 1 1 1 1 0 1 1 1 0 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 1 0 1 0 1 0 0 1
0 1 0 0 0 1 1 1 0 0 1 1 0 0 0 0 1 1 1 1 0 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

---

The minimum is $0$:

1 1 0 1 1 0 1 1 1 1 1 0 1 1 0 1 1 1 0 1 0 1 1 0 1 1 0 1 1 1 0 1 1 0 0 0 1 0 1 
0 0 1 0 1 0 1 1 1 1 1 1 1 1 0 1 1 0 0 0 1 0 1 1 0 0 1 1 0 1 1 0 1 1 1 1 1 1 1 
1 0 1 1 1 1 0 1 1 0 0 1 0 0 1 1 1 0 1 1 1 1 1 1 0 1 0 0 1 1 1 1 1 0 1 0 1 1 0 
1 0 1 0 1 0 1 0 1 0 0 1 0 0 1 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 
1 1 1 1 1 1 0 1 1 0 1 0 0 0 1 1 1 1 1 0 1 0 1 1 1 0 1 0 1 1 0 0 1 1 0 1 1 0 1 
1 1 1 0 0 1 1 1 0 1 0 1 1 1 1 1 0 0 1 1 1 0 0 1 0 1 0 1 1 0 1 0 1 1 1 0 1 1 1 
1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 0 1 0 1 1 0 0 1 0 0 1 1 1 1 1 1 0 0 0 1 1 0 1 0 
1 1 1 1 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 1 1 0 1 1 0 1 0 1 0 1 1 1 0 1 1 1 1 0 
1 0 0 0 1 1 0 1 1 1 0 1 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 0 
0 1 1 0 1 1 1 0 1 1 0 0 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 1 0 1 1 0 0 1 1 0 0 
0 1 1 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 0 1 1 1 1 1 0 0 0 0 1 1 1 0 1 1 0 0 1 1 
1 0 1 1 0 1 1 1 0 1 1 0 1 1 1 0 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 1 0 1 0 1 0 0 1 
0 1 0 1 0 1 1 1 0 0 1 1 0 0 0 1 1 1 1 1 0 1 0 0 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 

---

By focusing on the locations of the zeros, you can interpret a solution with objective value $0$ as a $(39,13,2)$-covering design with $13$ blocks, like the one found here: https://ljcr.dmgordon.org/show_cover.php?v=39&k=13&t=2