Optimal with IPOPTH, Infeasible with Gurobi

Hi, I’m working on a bilevel problem that is converted into an equivalent non-convex QCQP by applying strong duality to the lower-level problem (LP).

My code returns a (local) optimal solution with Ipopt solver (w/ ma86).
However, the same code returns infeasibility with Gurobi v9.0.2 (NonConvex=2).

This makes me wonder if the solution obtained with Ipopt is reliable.
Is there any possible reasoning for this? Any advice would be highly appreciated.

Thank you.

  [a076750e] CPLEX v0.6.5
  [336ed68f] CSV v0.6.2
  [aaaa29a8] Clustering v0.14.1
  [a93c6f00] DataFrames v0.21.3
  [31c24e10] Distributions v0.23.4
  [2e9cd046] Gurobi v0.8.1
  [f67ccb44] HDF5 v0.13.2
  [7073ff75] IJulia v1.21.2
  [b6b21f68] Ipopt v0.6.2
  [4076af6c] JuMP v0.21.3
  [e5e0dc1b] Juno v0.8.2
  [67920dd8] KNITRO v0.9.2
  [23992714] MAT v0.8.0
  [6405355b] Mosek v1.1.3
  [1ec41992] MosekTools v0.9.3
  [8bb1440f] DelimitedFiles

My experience has been that Ipopt is the most numerically stable non-linear solver. In this case, I would suspect there is a numerical issue in Guorbi. One useful diagnostic is to evaluate the feasibility of each constraint in the solution found by Ipopt.

3 Likes

@ccoffrin Thank you for your advice. Is there a smart way to check the feasibility?

The only way I know is to do something like: value.(g) - value.(d) == 0 for constraint @constraint(m, const1[t in 1:24], g[t] - d[t] == 0)
I just have too many constraints to do this manually.

Also, it is odd to me that Gurobi returns infeasibility rather than taking a longer time, as this means Gurobi fails to find any feasible point if I understand it correctly.

Also, it is odd to me that Gurobi returns infeasibility rather than taking a longer time, as this means Gurobi fails to find any feasible point if I understand it correctly.

I suggest you read https://www.gurobi.com/documentation/9.0/refman/num_grb_guidelines_for_num.html

Is there a smart way to check the feasibility?

Not really. You can use something like the following:

model = Model()
@variable(model, x >= 0)
@variable(model, y >= 0)
c = @constraint(model, x * y <= 1)
o = constraint_object(c)
value(o.func)  # left-hand side
o.set.upper  # right-hand side, if <=. Use .lower if >=, and .value if ==
1 Like

@odow Thank you so much for your comment!

1 Like