PolyJuMP makes a feasible model infeasible

[votroto]

When using PolyJuMP, the following model is reported as infeasible.

using JuMP, PolyJuMP, Gurobi

m = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))

@variable m -1 <= jump_c[1:5] <= 1
@variable m -1 <= jump_s[1:5] <= 1

c = jump_c
s = jump_s

# Solutions:
#c = [0.91165750302348, 0.9436635072926446, 0.575952449978958, 0.624599296124994, 0.9896625334970969]
#s = [0.4109508452126525, 0.3309066106987665, 0.8174831957681062, 0.7809454009597355, -0.14341572365716249]

@constraint m c .^ 2 .+ s .^ 2 .== 1
@constraint m (((c[1]*c[2] - 0.4608811486294164c[3]*s[4] - 0.5700328131603017c[5]*s[3] - 0.6801846504873382s[3]*s[5]) + (0.5700328131603017* ((c[3]*c[4]) * s[5]))) + (-0.6801846504873382* ((c[3]*c[4]) * c[5]))) == 0
@constraint m (((-0.7382491921224805c[3]*s[4] - 0.17909626144720603c[5]*s[3] + 0.6503173528871414s[3]*s[5] + s[2]) + (0.17909626144720603* ((c[3]*c[4]) * s[5]))) + (0.6503173528871414* ((c[3]*c[4]) * c[5]))) == 0
@constraint m (((0.5700328131603017c[3]*c[5] + 0.6801846504873382c[3]*s[5] - 0.4608811486294164s[3]*s[4] + s[1]) + (0.5700328131603017* ((c[4]*s[3]) * s[5]))) + (-0.6801846504873382* ((c[4]*c[5]) * s[3]))) == 0
@constraint m (((0.8018647772886565c[3]*c[5] - 0.3382841731078751c[3]*s[5] + 0.49252075810927465s[3]*s[4] - c[1]) + (0.8018647772886565* ((c[4]*s[3]) * s[5]))) + (0.3382841731078751* ((c[4]*c[5]) * s[3]))) == 0
@constraint m (((0.17909626144720603c[3]*c[5] - 0.6503173528871414c[3]*s[5] - 0.7382491921224805s[3]*s[4]) + (0.17909626144720603* ((c[4]*s[3]) * s[5]))) + (0.6503173528871414* ((c[4]*c[5]) * s[3]))) == 0
@constraint m 0.6503173528871414c[5]*s[4] + 0.17909626144720603s[4]*s[5] - c[2] + 0.7382491921224805c[4] == 0

optimize!(m)

but the commented lines for c and s contain the solution. If you un-comment them, the model becomes feasible. (You can convince yourself otherwise. The infeasibility of roughly 1e-16 is due to rounding)

Sorry this is the most I could clean it up for today, I'll try harder tomorrow.

[odow]

So I'm guessing the problem is related to

julia> using JuMP, PolyJuMP, Gurobi, Ipopt
       # Solutions:

julia> c_opt = [0.91165750302348, 0.9436635072926446, 0.575952449978958, 0.624599296124994, 0.9896625334970969]
5-element Vector{Float64}:
 0.91165750302348
 0.9436635072926446
 0.575952449978958
 0.624599296124994
 0.9896625334970969

julia> s_opt = [0.4109508452126525, 0.3309066106987665, 0.8174831957681062, 0.7809454009597355, -0.14341572365716249]
       # model = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))
5-element Vector{Float64}:
  0.4109508452126525
  0.3309066106987665
  0.8174831957681062
  0.7809454009597355
 -0.14341572365716249

julia> model = Model(Ipopt.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

julia> @variable(model, -1 <= c[i = 1:5] <= 1, start = c_opt[i])
5-element Vector{VariableRef}:
 c[1]
 c[2]
 c[3]
 c[4]
 c[5]

julia> @variable(model, -1 <= s[i = 1:5] <= 1, start = s_opt[i])
5-element Vector{VariableRef}:
 s[1]
 s[2]
 s[3]
 s[4]
 s[5]

julia> @constraints(model, begin
           c.^2 .+ s.^2 .== 1
           (((c[1]*c[2] - 0.4608811486294164c[3]*s[4] - 0.5700328131603017c[5]*s[3] - 0.6801846504873382s[3]*s[5]) + (0.5700328131603017* ((c[3]*c[4]) * s[5]))) + (-0.6801846504873382* ((c[3]*c[4]) * c[5]))) == 0
           (((-0.7382491921224805c[3]*s[4] - 0.17909626144720603c[5]*s[3] + 0.6503173528871414s[3]*s[5] + s[2]) + (0.17909626144720603* ((c[3]*c[4]) * s[5]))) + (0.6503173528871414* ((c[3]*c[4]) * c[5]))) == 0
           (((0.5700328131603017c[3]*c[5] + 0.6801846504873382c[3]*s[5] - 0.4608811486294164s[3]*s[4] + s[1]) + (0.5700328131603017* ((c[4]*s[3]) * s[5]))) + (-0.6801846504873382* ((c[4]*c[5]) * s[3]))) == 0
           (((0.8018647772886565c[3]*c[5] - 0.3382841731078751c[3]*s[5] + 0.49252075810927465s[3]*s[4] - c[1]) + (0.8018647772886565* ((c[4]*s[3]) * s[5]))) + (0.3382841731078751* ((c[4]*c[5]) * s[3]))) == 0
           (((0.17909626144720603c[3]*c[5] - 0.6503173528871414c[3]*s[5] - 0.7382491921224805s[3]*s[4]) + (0.17909626144720603* ((c[4]*s[3]) * s[5]))) + (0.6503173528871414* ((c[4]*c[5]) * s[3]))) == 0
           0.6503173528871414c[5]*s[4] + 0.17909626144720603s[4]*s[5] - c[2] + 0.7382491921224805c[4] == 0
       end);

julia> optimize!(model)
This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:       51
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       85

Exception of type: TOO_FEW_DOF in file "Interfaces/IpIpoptApplication.cpp" at line 662:
 Exception message: status != TOO_FEW_DEGREES_OF_FREEDOM evaluated false: Too few degrees of freedom (rethrown)!

EXIT: Problem has too few degrees of freedom.

This seems like a bug in the formulation that we are passing to Gurobi.

[votroto]

I would guess so. Hmm, for this smaller problem it is possible that the feasible set is zero-
dimensional, but this occurs also for problems (> 6-DOF IK) where that is known not to be the case.

Here is the reduced bug, as promised.

Removing the indicated parentheses makes the problem feasible.

using JuMP, PolyJuMP, Gurobi

m = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))

@variable m -1 <= c[3:5] <= 1
@variable m -1 <= s[3:5] <= 1

@constraint m s[4] == 0.7809454009597355
@constraint m 0.6796850426502927 - 0.4608811486294164 * s[3] * s[4] - 0.6731532644471365 * c[4] * s[3] == - 0.5700328131603017 * c[4] * s[3] * s[5]
@constraint m (-0.28365359534933055 - 0.1948355982905866 * s[5] + 0.49252075810927465 * s[3] * s[4]) + 0.8018647772886565 * c[4] * s[3] * s[5] == 0
#             ^                                                                                    ^
# Remove parentheses to make feasible

optimize!(m)

[odow]

Removing the parentheses is a very suspicious sign! Is it fixed by #112?

(-0.28365359534933055 - 0.1948355982905866 * s[5] + 0.49252075810927465 * s[3] * s[4]) makes a ScalarQuadraticFunction inside the ScalarNonlinearFunction, so I wonder if there's a similar problem.

[votroto]

I believe I tried it with the latest version directly from github main. I do not believe it is fixed.

[votroto]

Hmm... It looks as if the model is substituting the wrong variables.

If I model this:

m = Model(() -> PolyJuMP.QCQP.Optimizer(Gurobi.Optimizer()))

@variable m -1 <= v[1:4] <= 1

@constraint(m, v[3] == 0.780)
@constraint(m, 0.67 - 0.673*v[1]*v[2] + (0.57*v[1]*v[2]*v[4]) - 0.46*v[2]*v[3] == 0)
@constraint(m, -0.283 - 0.194*v[4] + 0.801*v[1]*v[2]*v[4] + 0.492*v[2]*v[3] == 0)

optimize!(m)

and inspect it with println(m.moi_backend.optimizer.model.model), depending on whether the parenthesized expr is on the lhs or rhs of the eq, the constraints change from

0.0 + 1.0 v[3] == 0.78
0.0 - 0.673 v[5] - 0.46 v[1]*v[3] + 0.57 v[2]*v[5] == -0.67
...

to

0.0 + 1.0 v[3] == 0.78
0.0 - 0.673 v[5] - 0.46 v[1]*v[2] + 0.57 v[3]*v[5] == -0.67
...                            ^

Based on the first constraint, the variable v[3] clearly refers to the same thing in all three models, and should therefore be associated with the coefficient 0.46, but in the second version, it is not.

I can't for the life of me figure out where the bug is. The code is a bit too complex for me.

[votroto]

BTW, there is a _to_polynomial! for ScalarQuadraticFunction in functions.jl and one with a union in nl_to_polynomial.jl (just the constraint on the dictionary changes). Are both necessary?

[odow]

Let me take a look.

[odow]

Here's a slightly simpler example:

julia> using JuMP, Gurobi, PolyJuMP

julia> function main()
           model = Model() do
               grb = Gurobi.Optimizer()
               MOI.set(grb, MOI.Silent(), true)
               PolyJuMP.QCQP.Optimizer(grb)
           end
           @variable(model, v[i = 1:3] >= i)
           @constraint(model, 12*v[1]*v[2] + (123*v[1]*v[2]*v[3]) == 1)
           @constraint(model, 3*v[3] + 123*v[1]*v[2]*v[3] == 2)
           optimize!(model)
           print(unsafe_backend(model).model)
       end
main (generic function with 1 method)

julia> main()
Feasibility

Subject to:

VariableIndex-in-GreaterThan{Float64}
 v[1] >= 1.0
 v[2] >= 2.0
 v[3] >= 3.0
 v[4] >= 3.0

ScalarQuadraticFunction{Float64}-in-EqualTo{Float64}
 0.0 + 1.0 v[4] - 1.0 v[1]*v[3] == 0.0
 0.0 + 12.0 v[4] + 123.0 v[2]*v[4] == 1.0
 0.0 + 3.0 v[2] + 123.0 v[2]*v[4] == 2.0

The affine terms in the bottom two constraints are wrong. It should be 12 * v[1] * v[2], but formulation is 12 * v[1] * v[3]and3.0 * v[2]instead of3.0 * v[3]`.

[votroto]

Is this right?

The right-hand-side seems un-ordered, that could cause issues.

I'm not sure why it's done this way. It seems that you just want the dictionary both ways: from MP to JuMP and from JuMP to MP. If I just pass the index_to_vars dict from final_touch to _add_variables! it's simple to get the inverse. It makes the original problem feasible, and your smaller problem seems to be correct too:

VariableIndex-in-GreaterThan{Float64}
 v[1] >= 1.0
 v[2] >= 2.0
 v[3] >= 3.0
 v[4] >= 2.0

ScalarQuadraticFunction{Float64}-in-EqualTo{Float64}
 0.0 + 1.0 v[4] - 1.0 v[1]*v[2] == 0.0
 0.0 + 12.0 v[4] + 123.0 v[3]*v[4] == 1.0
 0.0 + 3.0 v[3] + 123.0 v[3]*v[4] == 2.0

But I don't understand the code yet, so who knows... If it's just that, the functions could be simplified a bit.

[odow]

Want to make a PR?

[odow]

p.s., if you hadn't seen, I've been using one of your questions in talks: https://youtu.be/6q76umkG-34?si=KA-AzP8u_MbHF15u&t=215 😆

[votroto]

Yeah, sure, I'll open a PR.

Ooh, I feel honored :D

[votroto]

Hmm, this might take a while. It looks like there aren't all that many tests for the QCQP part. There are a few in a file called qcqp.jl but they don't run by default and some of them error due to missing MOI..

The error was in the _subs! function, where variables were getting reordered.

Also, the typing, seesh... starting with one(T), explicitly multiplying by floats, then converting back to T? I'm sure there is a reason, but can't see it now.

[votroto]

Actually, there is a comment

This will be refactored into a constraint bridge once jump-dev/MathOptInterface.jl#846 is done

the issue is closed.

[votroto]

While I'm at it, was this meant to check for divisibility, then use div_multiple, otherwise error with the "can't transform to QCQP" message?
Rational polynomials probably were not the intended outcome.

[blegat]

Yes, if we can avoid rational, it's indeed best. Otherwise there might be ways to handle rational polynomials but we can wait to have a use case before complicating the code for that