Why does OSQP return NaN for convex QP

Specific Domains Optimization (Mathematical)
,
1

I have a convex QP:

using JuMP, OSQP, Clarabel, Ipopt

A = [0.01; -0.01; 3.548942957597223e-5; -3.548942957597223e-5]
b = [6.998225528521202; 3.0017744714787984; 6.999996851250971; 3.1487490290693643e-6]
solver = optimizer_with_attributes(OSQP.Optimizer, "eps_abs" => 1e-4); #####Returns -3.702831462226849 
# solver = optimizer_with_attributes(OSQP.Optimizer, "eps_abs" => 1e-5); #####Returns NaN
# solver = optimizer_with_attributes(Ipopt.Optimizer);  #####Returns -0.08872357393992024
# solver = optimizer_with_attributes(Clarabel.Optimizer); #####Returns -0.08872357393992024

m = Model(solver)
set_silent(m)
@variable(m, w)
@expression(m, objExp, w^2 + 7.405663199206709 * w + 13.710961855021136)
@objective(m, Min, objExp)

@constraint(m, cstIn, A .* w .<= b)

optimize!(m)

OSQP returns either a wrong solution or NaN, depending on how I set the tolerances. Similar topics suggest nonconvexity of the QP but this is not the case here. So what is happening?

//Wasn’t sure if I should put it here or to OSQP github, let me know and I will move it.

2

For the "eps_abs" => 1e-4, OSQP is finding a solution that satisfies your tolerances (ish):

julia> optimize!(m)
-----------------------------------------------------------------
           OSQP v0.6.2  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 1, constraints m = 4
          nnz(P) + nnz(A) = 5
settings: linear system solver = qdldl,
          eps_abs = 1.0e-04, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, 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.3271e+01   1.52e-04   1.33e+00   1.00e-01   3.33e-05s
  25  -1.3711e+01   1.28e-04   2.50e-09   1.00e-01   4.44e-05s

status:               solved
number of iterations: 25
optimal objective:    -13.7110
run time:             4.93e-05s
optimal rho estimate: 4.62e+01


julia> primal_feasibility_report(m)
Dict{Any, Float64} with 1 entry:
  cstIn : -3.548942957597223e-5 w ≤ 3.1487490290693643e-6 => 0.000128263

For the NaN case, you must always check is_solved_and_feasible (or termination_status and primal_status) before querying the solution.

In your case, OSQP thinks the problem is infeasible:

julia> solution_summary(m)
* Solver : OSQP

* Status
  Result count       : 1
  Termination status : INFEASIBLE
  Message from the solver:
  "Primal_infeasible"

* Candidate solution (result #1)
  Primal status      : NO_SOLUTION
  Dual status        : INFEASIBILITY_CERTIFICATE
  Objective value    : 1.00000e+30
  Dual objective value : -3.14875e-06

* Work counters
  Solve time (sec)   : 5.22900e-05

Here’s a simpler example:

julia> using JuMP, OSQP

julia> begin
           model = Model(OSQP.Optimizer)
           set_attribute(model, "eps_abs", 1e-5)
           @variable(model, w)
           @objective(model, Min, w^2 + w)
           @constraint(model, 1e-2 * w >= -1)
           @constraint(model, 3e-5 * w >= 0)
           optimize!(model)
           solution_summary(model)
       end
-----------------------------------------------------------------
           OSQP v0.6.2  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 1, constraints m = 2
          nnz(P) + nnz(A) = 3
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, 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.9804e-01   2.18e-05   4.56e-01   1.00e-01   2.40e-05s
  25   1.0000e+30   1.50e-05   2.97e-09   1.00e-01   3.06e-05s

status:               primal infeasible
number of iterations: 25
run time:             3.45e-05s
optimal rho estimate: 8.48e+00

* Solver : OSQP

* Status
  Result count       : 1
  Termination status : INFEASIBLE
  Message from the solver:
  "Primal_infeasible"

* Candidate solution (result #1)
  Primal status      : NO_SOLUTION
  Dual status        : INFEASIBILITY_CERTIFICATE
  Objective value    : 1.00000e+30
  Dual objective value : -0.00000e+00

* Work counters
  Solve time (sec)   : 3.44610e-05

I don’t know if this is a “bug” in OSQP, or just the expected consequence of having constraints with small coefficients.

3

I’ve used OSQP (without the Jump wrapper) with few complaints for many problems. Feeding the OP’s problem directly into OSQP, I didn’t encounter the spurious feasibility issue, so I wonder how Jump is mangling the matrices.

In my experience, it is pretty clear that OSQP does not attempt to compensate for a poorly scaled problem (like the OP) whereas competitors seem to do so: this has implications for the interpretation of the tolerances. Here, scaling up the last (active) constraint yields a much more accurate solution.

1 Like