Presolve causes infeasibility error?

[devmotion]

When I tried to solve a trivial OT problem (10 sources with equal mass and a single target) with OptimalTransport.emd2 and Tulip, an error was thrown since the termination status was INFEASIBLE although the trivial solution is feasible and optimal. The following MWE seems to indicate that it is a problem with presolve: An optimal solution is found if presolve is disabled and otherwise the termination status is INFEASIBLE.

using Tulip
using JuMP

# Instantiate JuMP model
lp = Model(Tulip.Optimizer)

# Create variables
@variable(lp, x >= 0)
@variable(lp, y >= 0)
@variable(lp, z >= 0)

# Add constraints
@constraint(lp, row1, x + y + z == 1.0)
@constraint(lp, row2, x == 1/3)
@constraint(lp, row3, y == 1/3)
@constraint(lp, row4, z == 1/3)

# Set the objective
@objective(lp, Min, x + y + z)

# Set some parameters
# set_optimizer_attribute(lp, "Presolve_Level", 0)     # disable presolve

# Solve the problem
optimize!(lp)

# Check termination status
st = termination_status(lp)
println("Termination status: $st")

What is the recommended way for dealing with this issue? Should one disable presolve or are there some tolerances that should be adjusted? It seems a bit unfortunate if the default settings in OptimalTransport cause errors in these trivial but possibly numerically challenging cases, so I would be happy about any suggestions.

[mtanneau]

Thanks for reporting the issue. This is most likely a bug in the presolve.

Let me look into it.

[mtanneau]

OK, here is what happens:

• Presolve will fix x, then y, then z to 1/3, each time substracting 1/3 to the right-hand side of row1.
• This leaves row1 empty, but its RHS is now roughly 1e-16 because of round-off errors.
• There is no numerical tolerance when eliminating empty rows, so the problem is declared infeasible due to row1: 0 == 1e-16.

The good news is that it should not be too hard to fix 😅

Note that if you multiply all the rows by 3, then presolve returns the correct results

[devmotion]

Thanks for looking into this issue! I already assumed that it is at least partly caused by floating point issues. So if I understand you correctly, there's nothing that should be changed in our code but this will be addressed in the presolve algorithm at some point?

[mtanneau]

there's nothing that should be changed in our code but this will be addressed in the presolve algorithm at some point

That is correct. I "just" 😅 need to take an hour to add numerical tolerances into the presolve code.
I'll leave this issue open until then, and I can tag a release when I issue the fix.

[mtanneau]

@devmotion I've (finally!) added primal & dual tolerances in Presolve, which fixed the issue you were encountering.
I recommend you tag v0.7.5+

[devmotion]

Great, thanks a lot!