LQ problem infeasible OSQP after 0.6.1 update

[jd-lara]

after the update to 0.6.1 the following MWE returns infeasible.

gen =  [264.0
 331.0
 397.0
 462.0
 530.0]

price = [
9565.599999999999
 10791.66
 12036.52
 13340.8
 14748.7
]

m = Model(OSQP.Optimizer)
n_bp = length(price)
@variable(m, var_price[1:n_bp] >= 0)
@variable(m, quad_mult >= 0)
@variable(m, linear_mult >= 0)
@variable(m, intercept >= 0)
@constraint(m, [i in 1:n_bp], var_price[i] == intercept + linear_mult*gen[i] + quad_mult*gen[i]^2)
@objective(m, Min, sum((price[i] - var_price[i])^2 for i in 1:n_bp))
optimize!(m)

[jd-lara]

It seems to be related to the update to OSQP 0.6.2 binaries. With OSQP.jl v 0.6.0, I get this result.

[gbanjac]

It seems that the problem is close to being infeasible. If you decrease the primal infeasibility tolerance eps_prim_inf from the default value 1e-4 to 1e-6, OSQP should be able to solve the problem.

[RoyiAvital]

What's the reason the latest version is 0.6.1 and not 0.6.2?

[blegat]

What's your version of Julia and OSQP.jl ? Post Julia v1.3, this package should use OSQP_jll which is build using
https://github.com/JuliaPackaging/Yggdrasil/blob/master/O/OSQP/build_tarballs.jl
which seems to targer 0.6.2.
Pre Juila v1.3, you're using https://github.com/osqp/OSQP.jl/blob/master/deps/build.jl which seems to target 0.6.2 as well.
Of course, you need to be using the latest version of OSQP.jl

[RoyiAvital]

I looked at the release version of the repository. Not the binaries.

[odow]

If you decrease the primal infeasibility tolerance eps_prim_inf from the default value 1e-4 to 1e-6, OSQP should be able to solve the problem.

This isn't really true, because OSQP finds a solution that violates the constraints by quite a lot.

julia> using JuMP, OSQP

julia> function main()
           c = [264.0, 331.0, 397.0, 462.0, 530.0]
           p = [9565.599999999999, 10791.66, 12036.52, 13340.8, 14748.7]
           model = Model(OSQP.Optimizer)
           set_optimizer_attribute(model, "eps_prim_inf", 1e-6)
           @variable(model, x[1:5] >= 0)
           @variable(model, θ[1:3] >= 0)
           @constraint(model, [i=1:5], x[i] == θ[1] + c[i] * θ[2] + c[i]^2 * θ[3])
           @objective(model, Min, sum((p[i] - x[i])^2 for i = 1:5))
           optimize!(model)
           primal_feasibility_report(model)
       end
main (generic function with 1 method)

julia> main()
-----------------------------------------------------------------
           OSQP v0.6.2  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 8, constraints m = 13
          nnz(P) + nnz(A) = 33
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-06, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: off, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -1.3638e+06   3.73e-01   2.95e+04   1.00e-01   5.00e-05s
 200  -7.4832e+08   4.72e+02   2.96e+01   2.31e-05   1.09e-04s
 400  -7.4832e+08   3.20e+02   1.34e+01   2.13e-05   1.84e-04s
 600  -7.4830e+08   1.67e+01   1.43e+01   2.10e-04   2.56e-04s
 625  -7.4830e+08   1.24e+01   3.92e+00   2.10e-04   2.65e-04s

status:               solved
number of iterations: 625
optimal objective:    -748295532.8843
run time:             2.68e-04s
optimal rho estimate: 4.03e-04

Dict{Any, Float64} with 5 entries:
  x[3] - θ[1] - 397 θ[2] - 157609 θ[3] = 0.0 => 2.65857
  x[1] - θ[1] - 264 θ[2] - 69696 θ[3] = 0.0  => 1.63197
  x[4] - θ[1] - 462 θ[2] - 213444 θ[3] = 0.0 => 12.4135
  x[2] - θ[1] - 331 θ[2] - 109561 θ[3] = 0.0 => 1.76555
  x[5] - θ[1] - 530 θ[2] - 280900 θ[3] = 0.0 => 4.6714

I think this is just an upstream tolerance issue? The problem is pretty badly scaled.

julia> print(model)
Min x[1]² + x[2]² + x[3]² + x[4]² + x[5]² - 19131.199999999997 x[1] - 21583.32 x[2] - 24073.04 x[3] - 26681.6 x[4] - 29497.4 x[5] + 7.48339538956e8
Subject to
 x[1] - θ[1] - 264 θ[2] - 69696 θ[3] = 0.0
 x[2] - θ[1] - 331 θ[2] - 109561 θ[3] = 0.0
 x[3] - θ[1] - 397 θ[2] - 157609 θ[3] = 0.0
 x[4] - θ[1] - 462 θ[2] - 213444 θ[3] = 0.0
 x[5] - θ[1] - 530 θ[2] - 280900 θ[3] = 0.0
 x[1] ≥ 0.0
 x[2] ≥ 0.0
 x[3] ≥ 0.0
 x[4] ≥ 0.0
 x[5] ≥ 0.0
 θ[1] ≥ 0.0
 θ[2] ≥ 0.0
 θ[3] ≥ 0.0

I'd mark this as not a bug in OSQP.jl.

[imciner2]

It is probably related to osqp/osqp#346, which after a discussion between @gbanjac @bstellato and myself was traced back to a change in the tolerance checking code between 0.6.1 and 0.6.2. @gbanjac stated:

An infeasibility certificate y should satisfy A'y = 0 and supp(y) < 0. Previously, these conditions were implemented as ||A'y|| < eps and supp(y) < -eps. We then changed the second condition to supp(y) < eps as it makes more sense; the smaller the eps, the condition gets more similar to the true one. An issue is that for eps that is not that small, one can get false infeasibility detections.

This change makes the solver be slightly more pessimistic with problems that are close to infeasible instead of slightly optimistic. We talked about whether to revert that change, but couldn't really come to a final conclusion.