Presolve transforms your model into an equivalent model that has (in theory) the following properties:
- The presolved model is infeasible if and only if the original model is infeasible.
- The presolved model is unbounded if and only if the original model is unbounded.
- If the presolved model has an optimal solution, then its objective value is identical to the optimal objective value of the original model.
- Any feasible (or optimal) solution of the presolved model can be mapped to a feasible (or optimal) solution of the original model.
Moreover, the goal of presolve is that the presolved model will be smaller and easier to solve.
The optimize() function calls presolve to generate a presolved model p with the properties above and an "uncrush" mapping u: p -> m that can be used to map solutions of the presolved model p back to solutions of the original model m. Then, the presolved model is solved to find an optimal solution s' of p. Afterwards, this solution s' is mapped back to the original model by applying s := u(s').
Now, for numerical reasons, the properties above may not exactly be true, which means that s might be slightly infeasible or slightly sub-optimal for m. In this case (if s' came from a simplex solve or a barrier solve with crossover), we perform a few additional simplex iterations on the original model m starting from the basis that corresponds to s, in order to clean up infeasibilities and recover optimality.
So, what is the purpose of the Model.presolve() method in the Gurobi Python interface? It returns a model p that has the properties above, just like the one that is internally used by optimize(). It might be a little bit different from the one of optimize(), due to some technical reasons.
If you call p = m.presolve() followed by p.optimize(), you will get an optimal solution for p, and this solution will have the same objective value as an optimal solution for m. But because you do not have access to the mapping u: p -> m, you cannot convert this solution for p back into a solution for m. Many of the variables in p will be equivalent to the variables in m, which means that the mapping for these variables will be the identify function, but this is not true for all variables. For example, one presolve reduction is the "parallel columns" reduction. This merges two variables with identical objective and constraint coefficients into a single variable. For example, consider this problem:
max 3x1 + 3x2 + 4y + 5z
s.t. 2x1 + 2x2 + 3y <= 4
5x1 + 5x2 + 2z <= 7
2y + z <= 2
0 <= x1 <= 7
0 <= x2 <= 10
Variables x1 and x2 are "parallel". Presolve would replace these variables by a single variable x, which represents the sum of the two original variables:
max 3x + 4y + 5z
s.t. 2x + 3y <= 4
5x + 2z <= 7
2y + z <= 2
0 <= x <= 17
Now, let's assume we have a solution with x=15 for this model. Then, the uncrush function u could map this solution to x1=7, x2=8.
As a user, you would only see in this presolved model that there is some variable x that you weren't aware of and that this variable has a solution value of 15. You don't know how to get the solution values for your original variables x1 and x2.
Having said this, the presolve() method is mostly a debugging tool. Usually, you would write out this model to an *.lp file and inspect it manually.