Sorting using a MIP model

This is not, per se, very useful, but sorting a parameter inside a MIP model is not very easy for MIP solvers. Obviously, if you need a sorted parameter in the model, it is better to use a sorting algorithm. But useless models can still be interesting.

I use: 

    Input: a 1-dimensional parameter with values.
    Output: a 1-dimensional variable with the values sorted in ascending order.

We can implement this with a permutation matrix \(\color{darkred}X\), which is a permuted identity matrix. In a MIP context, this becomes a binary variable with some assignment constraints.

MIP Model for sorting \(p_i\)
\[\begin{align}\min\>&\color{darkred}z=0\\& \sum_i \color{darkred}x_{i,j} = 1&&\forall j\\ & \sum_j \color{darkred}x_{i,j} = 1&&\forall i\\ & \color{darkred}y_i = \sum_j \color{darkred}x_{i,j}\cdot\color{darkblue}p_j\\& \color{darkred}y_i \ge \color{darkred}y_{i-1}\\ & \color{darkred}x_{i,j} \in \{0,1\}\end{align}\]

If we want to sort \(n\) values, this model has \(n^2\) binary variables.

One possible improvement is to add the extra constraint \[\sum_i \color{darkred}y_i = \sum_i \color{darkblue}p_i\]

Let's compare a commercial solver (Cplex) with an open-source one (HiGHS). I use here \(n=50,100,200\). Model 1 is without the extra constraint while model 2 includes the extra constraint. The data \(\color{darkblue}p_i\) is random. A time limit of 1,000 seconds was used. Otherwise, the solvers used default settings.

----     96 PARAMETER results  solver performance

                        status        time       nodes  iterations

cplex.n=50 .model1     Optimal 2.970000E-1                    2408
cplex.n=50 .model2     Optimal 2.500000E-1                    2151
cplex.n=100.model1     Optimal           4                   14055
cplex.n=100.model2     Optimal           2                    7092
cplex.n=200.model1     Optimal         227        3043     1684124
cplex.n=200.model2     Optimal          18                   10155
highs.n=50 .model1     Optimal           2
highs.n=50 .model2     Optimal           2
highs.n=100.model1   TimeLimit        1001        9587     4955336
highs.n=100.model2   TimeLimit        1000       15221     5373228

Some conclusions:

• Cplex is much better on the larger models
• HiGHS is very promising on the \(n=50\) model: it solves this in the presolver (zero nodes, zero LP iterations). But it runs out of steam for large models. I would expect that if the presolver can solve the problem for \(n=50\), it would also be able to do this for \(n=100\). As usual, my predictions are not that good.
• The extra constraint is useful for Cplex, especially for the larger models. Both in terms of time and LP iterations.

The Xpress solver was complaining about numerical issues:

Numerical issues encountered:

   Dual failures    :    831 out of      5345 (ratio: 0.1555)

   Singular bases   :      1 out of     50821 (ratio: 0.0000)

I am a little bit surprised by this. Scaling is the most notorious reason for numerical problems in LP/MIP models. Since this model does not have big-M constraints, scaling is excellent. Of course, I don't how the cuts added by the solver look like. I don't really know what to do about this. Maybe the absence of an objective is a problem.

Here, we sorted a parameter. That is relatively easy. Sorting a variable inside a model is more of a challenge. The constraint \[\color{darkred}y_i = \sum_j \color{darkred}x_{i,j}\cdot\color{darkred}p_j\] is now quadratic. The products \(\color{darkred}x_{i,j}\cdot\color{darkred}p_j\)  (a binary times a continuous variable) can be linearized, but that will give us some big-M constraints. Indicator constraints are an option if no good bounds on \(\color{darkred}p_j\) are available. This is again a case where using quadratic formulations can make the modeling way easier.