Issue with OSQP Returning Negative Values Despite Non-Negative Constraints

[baharehhj]

Hi everyone,

I hope you're all doing well. I'm reaching out because I've encountered an issue while working with the OSQP solver for a quadratic programming problem.

In brief, my problem involves minimizing a quadratic objective function subject to linear inequality constraints. To address this, I've integrated the OSQP solver with the Eigen library and OSQP-Eigen in C++.

Here's a quick overview of the problem components:

Objective Function: Quadratic, defined by matrix H and vector f.
Inequality Constraints: Defined by the sparse matrix Ac and vector dc.
Bounds: Lower bounds lb and upper bounds ub.
Non-Negativity Constraints: I've implemented constraints to ensure non-negativity for certain variables. To enforce this, I've constructed a diagonal matrix (identity matrix) and concatenated it with the inequality constraints matrix Ac.
However, after setting up the problem and solving it using OSQP, I'm encountering negative values in the solution for some variables. This is unexpected, as the problem explicitly requires non-negative solutions.

I've reviewed the problem setup, adjusted solver parameters, and examined the solution post-solve, but I haven't been able to identify the root cause of the negative values.

For those interested in diving into the code, I've shared the relevant snippets in the GitHub repository. Here's the link to the specific file.

Thank you all in advance for any insights or assistance you can provide.

[goulart-paul]

How negative are the entries in the solution? If you are getting negative values on the order of -1e-5 or so (depending on solver parameters), then this is expected behaviour since the solver only produces solutions to within the tolerances it has set.

If you are getting values that are so negative that they can't be explained by solver precision, then it's likely that either 1) there is a error in the problem setup or 2) the solver is reporting that the problem is infeasible, in which case the solution vector it returns has a different interpretation.

[baharehhj]

Thank you for your quick answer.

The smallest negative number in the solution is -0.01.

With these settings:

solver.settings()->setRho(0.1);
solver.settings()->setMaxIteration(8000);
solver.settings()->setWarmStart(true);
solver.settings()->setVerbosity(true);
solver.settings()->setAdaptiveRho(true);
solver.settings()->setLinearSystemSolver(0);
solver.settings()->setScaling(2);
solver.settings()->setAbsoluteTolerance(1.0e-5);
solver.settings()->setRelativeTolerance(1.0e-5);

[goulart-paul]

It is still hard to give a clear answer here I think. The solver terminates on a "relative accuracy" check, so whether or not it considers "-0.1" to be a large violation depends on what the norm of the solution x is in the first place. I would suspect that either that norm is large, or the solver is actually reporting that the solution is not feasible. Did you check both of those conditions?

[baharehhj]

Thank you for your continued assistance
The problem has been solved, and it is feasible. A part of the solution, as I wrote earlier, is as follows:

problem:  variables n = 721, constraints m = 2370
          nnz(P) + nnz(A) = 470615
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-05,
          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 = 8000
          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  -3.9093e+02   7.79e+00   6.40e-01   1.00e-01   1.83e-01s
 200  -6.6152e+04   3.33e+00   1.01e-03   2.54e-05   9.11e-01s
 400  -6.6163e+04   2.17e-01   3.89e-04   2.54e-05   1.30e+00s
 600  -6.6168e+04   1.18e+00   1.26e-04   4.25e-06   1.79e+00s
 800  -6.6170e+04   1.97e-01   6.93e-05   4.25e-06   2.15e+00s
1000  -6.6171e+04   8.81e-02   5.10e-05   4.25e-06   2.57e+00s
1200  -6.6171e+04   5.90e-02   4.46e-05   4.25e-06   2.94e+00s
1400  -6.6171e+04   4.23e-02   3.67e-05   4.25e-06   3.30e+00s
1600  -6.6171e+04   3.10e-02   2.87e-05   4.25e-06   3.70e+00s
1800  -6.6171e+04   2.25e-02   1.89e-05   4.25e-06   4.07e+00s
1925  -6.6171e+04   1.96e-02   1.48e-05   4.25e-06   4.29e+00s

status:               solved
number of iterations: 1925
optimal objective:    -66171.4063
run time:             4.29e+00s
optimal rho estimate: 1.07e-06

solution:
-1.66835e-07
-1.35141e-07
 2.29064e-06
-2.13231e-06
 4.49156e-07
 -6.5708e-07
-5.01667e-07
-4.49771e-07
-3.01505e-05
 3.97286e-05
 3.16376e-05
     31.6182
     71.5931
     42.5751
     9.96366
-3.84658e-07
 -4.6553e-07
     1547.52
     1428.15
     2658.23
     37.0342
 -9.1152e-05
     227.749
     77.6382
-4.56218e-08
      1.1573
-4.83679e-07
  0.00061445
     4832.22
     2667.62
     3007.15
     1339.03
     1948.28
     300.566
-0.000457403
-4.04041e-05
     3.04869
-4.75612e-07
     1587.39
      2943.4
      405.77
     4168.62
     1967.75
     2560.26
     1523.52
     1882.96
     159.161
     12.5638
     1.68854
     10.6766
 5.87291e-06
     2679.55
     2157.07
     498.351
     2555.83
     1589.89
     1954.24
     1111.21
     1651.56
     467.314
     34.3023
  0.00129917
      1.9675
 1.06878e-05
     2335.01
      717.91
  -0.0052243
   -0.00116751
-0.000273156
 1.67115e-05
     49.5698
     42.7964
      1788.7
     1614.46
  -0.0131135
     2109.85
     99.6026
     1227.37
  0.00232391
 0.000872348

[goulart-paul]

The convergence criteria for the solver is in section 3.4 of the paper here. The problem you are having is that the convergence check compares the error residuals to the infinity norm of the product $Ax$. It seems likely (based on the solution you show above) that $|Ax| \ge {}$ 1e+3. The error residuals themselves (thinking only about the nonnegativity constraint) look to have norm of about $|r| ={}$ 1e-3 (this is the size of the largest magnitude negative number).

The solver is therefore satisfying the convergence test in (26) of the paper cited above, i.e. $|r| \le \epsilon_{abs} + \epsilon_{rel} \max({|Ax|,|z|} )$, since you have $\epsilon_{rel}$ set to 1e-5.

If you want higher accuracy, you will need to tighten $\epsilon_{rel}$. You could also try setting polish = true, which will give a very high accuracy solution if it works but is not guaranteed to do so.

[baharehhj]

When I solve this optimization problem using the MATLAB function quadprog, I obtain a solution where none of the numbers are negative. In my problem, it is necessary to have non-negative numbers in the solution. Therefore, I have to change any negative values to zero. Does this alteration make a significant difference?

When I set the relative tolerance and absolute tolerance to 1e-3, I observe fewer negative numbers, for example, 17. However, when I set them to 1e-5, I encounter 100 negative numbers. Furthermore, when I set the tolerances to 1e-7, I still get about 100 negative numbers, with the smallest negative number being -0.0001.

If I set polish to True, I get the following result:
Solution polish is unsuccessful.

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

iter  objective    pri res    dua res    rho        time
   1  -9.0300e+01   1.68e+00   5.18e-01   5.00e-01   1.95e-01s
 200  -6.6136e+04   1.26e+00   3.50e-03   4.14e-04   7.98e-01s
 400  -6.6167e+04   2.40e-01   8.51e-04   6.93e-05   1.30e+00s
 600  -6.6175e+04   1.05e-01   5.45e-04   6.93e-05   1.67e+00s
 800  -6.6179e+04   7.00e-01   3.11e-04   1.28e-05   2.18e+00s
1000  -6.6182e+04   2.20e-01   1.78e-04   1.28e-05   2.56e+00s
1200  -6.6184e+04   1.25e-01   1.24e-04   1.28e-05   2.94e+00s
1400  -6.6185e+04   8.45e-02   1.01e-04   1.28e-05   3.31e+00s
1600  -6.6185e+04   5.64e-02   8.47e-05   1.28e-05   3.69e+00s
1800  -6.6185e+04   4.01e-02   7.10e-05   1.28e-05   4.08e+00s
2000  -6.6186e+04   3.01e-02   5.96e-05   1.28e-05   4.46e+00s
2200  -6.6186e+04   2.29e-02   5.15e-05   1.28e-05   4.83e+00s
2400  -6.6186e+04   1.47e-02   4.49e-05   1.28e-05   5.22e+00s
2600  -6.6186e+04   1.22e-02   3.94e-05   1.28e-05   5.60e+00s
2800  -6.6186e+04   2.10e-01   2.35e-05   2.44e-06   6.09e+00s
3000  -6.6186e+04   1.34e-01   1.45e-05   2.44e-06   6.47e+00s
3200  -6.6186e+04   8.75e-02   8.98e-06   2.44e-06   6.84e+00s
3400  -6.6186e+04   6.00e-02   5.57e-06   2.44e-06   7.22e+00s
3600  -6.6186e+04   4.43e-02   4.13e-06   2.44e-06   7.60e+00s
3800  -6.6186e+04   3.24e-02   3.06e-06   2.44e-06   7.98e+00s
4000  -6.6186e+04   2.47e-02   2.16e-06   2.44e-06   8.36e+00s
4200  -6.6186e+04   1.94e-02   1.51e-06   2.44e-06   8.74e+00s
4400  -6.6186e+04   1.57e-02   1.06e-06   2.44e-06   9.11e+00s
4600  -6.6186e+04   1.29e-02   7.42e-07   2.44e-06   9.49e+00s
4800  -6.6186e+04   1.14e-02   5.21e-07   2.44e-06   9.87e+00s
5000  -6.6186e+04   1.04e-02   3.67e-07   2.44e-06   1.03e+01s
5200  -6.6186e+04   9.44e-03   2.61e-07   2.44e-06   1.06e+01s
5400  -6.6186e+04   8.52e-03   1.86e-07   2.44e-06   1.10e+01s
5550  -6.6186e+04   7.88e-03   1.44e-07   2.44e-06   1.13e+01s

status:               solved
solution polish:      unsuccessful
number of iterations: 5550
optimal objective:    -66186.2536
run time:             1.14e+01s
optimal rho estimate: 6.29e-06

[goulart-paul]

There's no problem with rounding the values up to zero yourself. The issue is that the solver can't do that for you because it treats your constraint in a general $l \le Ax \le u$ form. It's obvious how to round if $A = I$, but not obvious at all in the general case.

If you have a collection of constraints -- $Ax \le b$, $x \ge 0$, say -- then you can round up the $x$ values to satisfy the second set of conditions. You may end up violating the first one though. Or, in the solver's terms, you mean end up with a violation residual greater than what the solver was aiming for.

Note that the number of constraint violations is not a metric used by the solver. It only uses the relative size of the residuals. So when you say:

Furthermore, when I set the tolerances to 1e-7, I still get about 100 negative numbers, with the smallest negative number being -0.0001.

that means that the solver has produced a valid solution according to its convergence criteria, assuming that the norm of $x$ is as in the earlier discussion.

Polishing would likely have resolved the issue for you if successful, but it doesn't always work. When it doesn't work, we just return the solution we computed prior to the polishing attempt.

[baharehhj]

Thank you