MILP formulation and optimization

Question

For $i=1, \dotsc, K$, we have $n_i$ ordered real numbers:

$$ x_i(1) \leq x_i(2) \leq \dotsc \leq x_i(n_i) $$

I want to solve the following optimization problem:

\begin{align} \mathrm{maximize} \; \sum_{i=1}^Kr_i \;\; s.t. \\ r_i \in \{1,\dotsc, n_i\}\\ \sum_{i=1}^Kr_i \leq c\sum_{i=1}^K x_i(r_i) \end{align}

Now, I had the idea of expressing this problem as a MILP as follows:

For each $i$, define binary variables which satisfy the linear constraints:

\begin{align} z_i(1) \geq z_i(2) \geq \dotsc \geq z_i(n_i) \in \{0,1\} \end{align}

Also define $y_i(k) = x_i(k)-x_i(k-1), k = 1,\dotsc, n_i$, where $x_i(0)=0$.

Then we can solve the MILP:

\begin{align} \mathrm{maximize} \; \sum_{i=1}^K \sum_{j=1}^{n_i} z_i(j) \;\; s.t. \\ z_i(1) \geq z_i(2) \geq \dotsc \geq z_i(n_i) \in \{0,1\} \\ \sum_{i=1}^K \sum_{j=1}^{n_i} z_i(j) \leq c\sum_{i=1}^K \sum_{j=1}^{n_i} y_i(j)z_i(j) \end{align}

Now I have the following questions:

1. Is there any obvious easier way of modelling this problem? e.g. as a MILP with less integer variables or using another optimization methodology.

2. Is there any simple way to intervene with the solvers in order to improve the solution speed of this problem?

3. In most of my problems, the number of binary variables is between 20000-100000. For this reason, I was actually very surprised that using the Symphony solver, those problems get solved in less than a minute! Should these problem sizes not be prohibitive?

4. Based on prior information, I also tried to remove the highest $x_i(j)$ values, thus reducing the total number of binary variables. E.g. reducing the number of variables by 50% often lead to substantial speed gains. On the other hand, it was often the case that when I reduced the number of variables to 10%, the solver would actually take longer than with all the variables (even though it would eventually return the same solution). How can this be explained?

Answer 1

1. Is there any obvious easier way of modelling this problem? e.g. as a MILP with less integer variables or using another optimization methodology.

The first tactic that immediately comes to mind is, for each $i$, modeling the set $\{z_{i,1}, \ldots, z_{i,n_{i}}\}$ as a SOS1 set. Your current formulation essentially does the same thing.

2. Is there any simple way to intervene with the solvers in order to improve the solution speed of this problem?

Aside from reformulating the problem to use SOS1 sets, you could use a better solver. For instance, you could use either Gurobi or CPLEX, both of which regularly compete for best performing solver over test beds of MILP and LP instances. Both of these solvers are roughly 10x faster than any open source MILP solver. Changing your solver is probably the fastest way to achieve some sort of speed-up.

In addition, you could look to see if you could reformulate your problem so that the LP relaxations are tighter, or introduce more and better cuts that exclude nonintegral regions of your feasible region. A senior colleague of mine has pointed out that adding custom cuts in CPLEX will disable its very robust presolve routines, along with its cut heuristics; I have not seen this behavior in practice, but the mere possibility is worth keeping an eye on if you choose to go this route.

You might also be able to come up with a polynomial time approximation scheme that is either "good enough", or maybe it even generates feasible solutions for your formulation. These feasible solutions would provide upper bounds on the optimal objective function value, which would then be useful in fathoming nodes in a branch-and-bound tree. Your formulation looks vaguely like a knapsack problem, but it doesn't have the right structure to be a knapsack problem. Nevertheless, you might be able to impose some structure like that to obtain a decent polynomial time approximation algorithm (or even a solution algorithm).

3. In most of my problems, the number of binary variables is between 20000-100000. For this reason, I was actually very surprised that using the Symphony solver, those problems get solved in less than a minute! Should these problem sizes not be prohibitive?

Structure is everything. Since each group of $z_{i,1}, \ldots, z_{i,n}$ is a SOS1 set, a branch-and-bound implementation only needs to examine $\prod_{i = 1}^{K}n_{i}$ possible combinations of binary variables instead of $2^{\sum_{i=1}^{K}n_{i}}$ variables. The former expression grows more slowly than the latter. I suspect CPLEX or Gurobi would preprocess this problem, deduce this structure, and employ effective cuts to solve this problem quickly (relative, to say, a random MILP instance with 20000-100000 binary variables and no continuous variables).

Answer 2

I will try to give one quick answer to my 2. question, based on Geoff's accepted answer. Geoff suggested to use a polynomial approximation scheme to get feasible solutions (which in turn can speed-up the branch-and-bound process). One such scheme which will produce a quite high quality solution in $O(nlogn)$ operations ($n$ being the number of binary variables in the original problem) is the following:

1. Merge the groups, i.e. discard all the group information and pool all $x_i(j)$'s together.

2. Sort the above list.

3. Now the problem has been reduced to the case $K=1$. Just one run through the list is necessary in order to determine the best solution of this problem which has objective value $r$.

4. We now look at the $r$-lowest $x_i(j)$'s. By counting how many of these correspond to each of the initial groups, we can get the $r_i$ values. This is a feasible solution and in most cases of quite high quality.