Gurobi sometimes fails to return duals in a QCP

[odow]

@adow031 has seen an example where Gurobi solves a QCP, experiences numerical error and fails to return the duals, but DualStatus is FEASIBLE_POINT.

The error is

Barrier performed 54 iterations in 0.08 seconds
Sub-optimal termination - objective 2.48377287e+05

Warning: failed to compute QCP dual solution due to inaccurate barrier solution
         Try decreasing BarQCPConvTol for more accuracy

(This issue is to remind me that the problem exists. I'll try to find a reproducible example that isn't "run this large SDDP problem". @adow031 will post a .mof.json file when he next triggers it.)

I guess it isn't sufficient to just consider the SolCount. We should probably check that we can query RC.

Gurobi.jl/src/MOI_wrapper/MOI_wrapper.jl

Lines 2833 to 2836 in 32b15ed

	# Check SolCount explicitly instead of ResultCount to avoid returning 1 when
	# there is a certiticate.
	valueP = Ref{Cint}()
	ret = GRBgetintattr(model, "SolCount", valueP)

[odow]

@haoxiangyang89 found the likely culprit in odow/SDDP.jl#477: the dual isn't defined at the point of a SecondOrderCone. I'll see if I can reproduce.

@adow031 / @DrAndyPhilpott you should add a constraint to avoid the case 0 >= 0 in a quadratic constraint.

[adow031]

@haoxiangyang89 found the likely culprit in odow/SDDP.jl#477: the dual isn't defined at the point of a SecondOrderCone. I'll see if I can reproduce.

@adow031 / @DrAndyPhilpott you should add a constraint to avoid the case 0 >= 0 in a quadratic constraint.

We found that running our model with CPLEX didn't have any errors; so it may be a Gurobi specific issue. However, running SDDP.jl with a QP took a very long time, so we have gone back to using a piecewise-linear approximation.

[MathiasPBerger]

Hi guys,

I don't know if this is still relevant but I just ran into the same problem on a fairly simple QCQP example, please see here.

The example was run under Julia 1.6.5 and uses JuMP v0.23.2 and Gurobi v0.11.0. Termination status is 'OPTIMAL' and dual status is 'FEASIBLE POINT' but the 'Pi' attribute cannot be retrieved.

[odow]

Oh great! I've been struggling to reproduce this.

Can you try adding constraints such that

(p_max[g] - p[g] + sum(alpha[g, i]*mu[i] for i = 1:n_w)) >= 0.001

My hunch is that the problem occurs at the point of the second-order cone.

[MathiasPBerger]

Will do this and let you know.

[MathiasPBerger]

The error indeed disappears when this constraint is added for each generator (g index). The optimal objectives are almost equal in both cases (with and without the constraint that you proposed).

I'm a bit perplexed by the fact that the matrix I'm using to build the second-order cone constraint is positive definite and the vector that multiplies it has nonzero entries in the solution (including in the case where the error occurs), so I don't exactly see how we could end up with a 0 >= 0 situation.

[odow]

If you don't include the constraints, what is value.(max_prod) and value.(min_prod) of the optimal solution?

[MathiasPBerger]

This is what I get:

max_prod: [50.27933145432515, 1.0638298606949077, 1.0638298606977323]
[68.25933392790412, 0.7978734083340456, 0.7978734083360246]
[79.04733539688482, 0.638296730977071, 0.6382967309721246]

min_prod: [71.31205574336671, 1.063829860695034, 1.0638298606982575]
[53.33205327073633, 0.7978734083327054, 0.7978734083346006]
[42.54405178742139, 0.6382967309756538, 0.6382967309741281]

I solved the same problem with Mosek and no error occurred, which supports the idea that this is a Gurobi-related issue.

[odow]

This is what I get

Interesting. We obviously need to dig a bit deeper to understand what's going on. I don't really have time to do this, unfortunately. cc @simonbowly

I solved the same problem with Mosek and no error occurred, which supports the idea that this is Gurobi-related issue

The problem is that Gurobi doesn't actually solve a second-order cone, it solves t^2 >= x^2 + y^2.

[MathiasPBerger]

Hi guys,

In case somebody picks it up later, I updated the formulation a bit and turned the original SOC constraints into quadratic constraints of the form x^T A x <= y^2 (along with a couple of affine equality constraints), where y is a non-negative scalar variable and A is positive definite.

Gurobi manages to solve the updated formulation without any problem, including for parameter values that proved problematic for the original SOCP formulation. The objectives and primal solutions are the same for both formulations.

Best,

Mathias