Unable to solve LP accurately

I’ve just started using OSQP and I have a small LP test problem that seems to cause it trouble and I’m trying to figure out why. Any insight is appreciated. I’m using the MATLAB interface on macOS.

Using the default options, it claims it solves successfully (in 3425 iterations), but the solution is significantly different from the consistent solution found by multiple other solvers. Attempting to tighten the termination tolerances or turn on polish results in failure.

The example problem, and results from OSQP and quadprog() are below.

Code for example problem:

% problem data
P = sparse(16, 16);
q = [zeros(11,1); 2403.5; 1; 1; 1000; 1];
i = [1 1 1 2 2 2 3 3 3 4 4 4 4 5 5 5 6 6 6 6 7 7 7 8 8 8 8 9 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 21 22 22 22 23 24 25 26 27 28]';
j = [1 4 10 2 8 12 3 6 11 1 4 5 9 4 5 6 3 5 6 7 6 7 8 2 7 8 9 4 8 9 5 6 6 7 7 8 2 8 8 9 4 9 10 13 10 13 10 13 11 14 11 14 10 12 15 11 12 16 1 10 11 12 15 16]';
a = [17.3611111 -17.3611111 -1 16 -16 -1 17.0648464 -17.0648464 -1 -17.3611111 39.9953822 -10.8695652 -11.7647059 -10.8695652 16.7519182 -5.88235294 -17.0648464 -5.88235294 32.8678343 -9.92063492 -9.92063492 23.8095238 -13.8888889 -16 -13.8888889 36.100069 -6.21118012 -11.7647059 -6.21118012 17.975886 5.88235294 -5.88235294 9.92063492 -9.92063492 13.8888889 -13.8888889 -16 16 6.21118012 -6.21118012 -11.7647059 11.7647059 2500 -1 3000 -1 3500 -1 1500 -1 2000 -1 -1 1 -1 -1 1 -1 1 1 1 1 1 1]';
A = sparse(i, j, a, 28, 16);
l = [0 0 0 0 -0.9 0 -1 0 -1.25 -1.5 -0.4 -2.5 -2.5 -2.5 -2.5 -Inf -Inf -Inf -Inf -Inf -Inf -Inf 0 0.9 0.1 0.1 0 0]';
u = [0 0 0 0 -0.9 0 -1 0 -1.25 1.5 0.4 2.5 2.5 2.5 2.5 0 500 1500 0 1000 0 0 0 2.5 2.7 3 Inf Inf]';

% solve via OSQP
prob = osqp;
prob.setup(P, q, A, l, u, opt);
res = prob.solve();
lb_res = A * res.x - l;
ub_res = u - A * res.x;

% solve via quadprog
[x, f, eflag] = quadprog(P, q, [A;-A], [u;-l]);

% compare results
fprintf('max constraint violation = %g\n', abs(min([lb_res; ub_res])))
fprintf('max difference in x = %g\n', norm(x - res.x, Inf));
fprintf('obj val difference = %g\n', norm(f - res.info.obj_val, Inf));

Results:

-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 16, constraints m = 28
          nnz(P) + nnz(A) = 64
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, 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.1654e+04   5.87e+00   7.48e+05   1.00e-01   5.58e-05s
 200   6.1629e+03   5.35e-01   1.54e+03   4.92e-02   1.68e-04s
 400   7.8530e+03   3.17e-01   4.19e+02   4.92e-02   2.72e-04s
 600   7.1161e+03   8.19e-02   6.35e+02   4.92e-02   3.76e-04s
 800   7.2214e+03   3.16e-01   8.26e+01   5.52e-03   4.82e-04s
1000   6.8787e+03   1.76e-01   7.77e+01   5.52e-03   5.88e-04s
1200   6.7038e+03   1.87e+00   8.86e+01   5.52e-03   6.93e-04s
1400   6.5888e+03   2.12e-01   4.53e+01   5.52e-03   7.97e-04s
1600   6.9733e+03   2.01e-02   7.19e+01   5.52e-03   9.01e-04s
1800   7.2828e+03   1.32e-01   5.35e+01   5.52e-03   1.00e-03s
2000   7.0420e+03   5.10e-02   4.84e+01   5.52e-03   1.11e-03s
2200   6.8012e+03   7.67e-02   1.34e+01   5.52e-03   1.21e-03s
2400   6.9411e+03   4.85e-02   5.56e+01   5.52e-03   1.32e-03s
2600   7.2646e+03   1.14e-01   3.67e+01   5.52e-03   1.42e-03s
2800   7.0907e+03   4.50e-02   5.54e+01   5.52e-03   1.52e-03s
3000   6.7481e+03   1.40e-01   1.47e+01   5.52e-03   1.63e-03s
3200   6.9816e+03   3.59e-02   4.20e+01   5.52e-03   1.73e-03s
3400   7.2479e+03   9.89e-02   3.30e+01   5.52e-03   1.84e-03s
3425   7.2426e+03   8.85e-02   1.51e+00   5.52e-03   1.86e-03s

status:               solved
number of iterations: 3425
optimal objective:    7242.6069
run time:             1.86e-03s
optimal rho estimate: 3.56e-02

Warning: Quadratic matrix (H) is empty or all zero. Your problem is linear.
Consider calling LINPROG instead. 
> In quadprog (line 336)
  In osqp_lp_fail (line 14) 

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

max constraint violation = 0.0446026
max difference in x = 421.842
obj val difference = 225.373

I couldn’t get your code to run exactly as written, because opt is not declared anywhere. Perhaps it is declared as an empty array somewhere, and that is the cause of the problems?

I tidied it up and then ran with polishing enabled, which seemed to give exactly the same solution as quadprog. The code as I ran it (some cosmetic edits relative to your example):

% problem data

P = sparse(16, 16);

q = [zeros(11,1); 2403.5; 1; 1; 1000; 1];

i = [1 1 1 2 2 2 3 3 3 4 4 4 4 5 5 5 6 6 6 6 7 7 7 8 8 8 8 9 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 21 22 22 22 23 24 25 26 27 28]';

j = [1 4 10 2 8 12 3 6 11 1 4 5 9 4 5 6 3 5 6 7 6 7 8 2 7 8 9 4 8 9 5 6 6 7 7 8 2 8 8 9 4 9 10 13 10 13 10 13 11 14 11 14 10 12 15 11 12 16 1 10 11 12 15 16]';

a = [17.3611111 -17.3611111 -1 16 -16 -1 17.0648464 -17.0648464 -1 -17.3611111 39.9953822 -10.8695652 -11.7647059 -10.8695652 16.7519182 -5.88235294 -17.0648464 -5.88235294 32.8678343 -9.92063492 -9.92063492 23.8095238 -13.8888889 -16 -13.8888889 36.100069 -6.21118012 -11.7647059 -6.21118012 17.975886 5.88235294 -5.88235294 9.92063492 -9.92063492 13.8888889 -13.8888889 -16 16 6.21118012 -6.21118012 -11.7647059 11.7647059 2500 -1 3000 -1 3500 -1 1500 -1 2000 -1 -1 1 -1 -1 1 -1 1 1 1 1 1 1]';

A = sparse(i, j, a, 28, 16);

l = [0 0 0 0 -0.9 0 -1 0 -1.25 -1.5 -0.4 -2.5 -2.5 -2.5 -2.5 -Inf -Inf -Inf -Inf -Inf -Inf -Inf 0 0.9 0.1 0.1 0 0]';

u = [0 0 0 0 -0.9 0 -1 0 -1.25 1.5 0.4 2.5 2.5 2.5 2.5 0 500 1500 0 1000 0 0 0 2.5 2.7 3 Inf Inf]';

 

% solve via OSQP

prob = osqp;

prob.setup(P, q, A, l, u,'polish',1);

results = prob.solve();

 

lb_res = A * results.x - l;

ub_res = u - A * results.x;

 

% solve via quadprog

[x, f, eflag] = quadprog(P, q, [A;-A], [u;-l]);

 

% compare results

fprintf('max constraint violation = %g\n', abs(min([lb_res; ub_res])))

fprintf('max difference in x = %g\n', norm(x - results.x, Inf));

fprintf('obj val difference = %g\n', norm(f - results.info.obj_val, Inf));

And the result:

-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 16, constraints m = 28
          nnz(P) + nnz(A) = 64
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, 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: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -1.1654e+04   5.87e+00   7.48e+05   1.00e-01   8.99e-05s
 200   6.1629e+03   5.35e-01   1.54e+03   4.92e-02   2.36e-04s
 400   7.8530e+03   3.17e-01   4.19e+02   4.92e-02   3.70e-04s
 600   7.1161e+03   8.19e-02   6.35e+02   4.92e-02   5.02e-04s
 800   7.2214e+03   3.16e-01   8.26e+01   5.52e-03   6.33e-04s
1000   6.8787e+03   1.76e-01   7.77e+01   5.52e-03   7.63e-04s
1200   6.7038e+03   1.87e+00   8.86e+01   5.52e-03   8.93e-04s
1400   6.5888e+03   2.12e-01   4.53e+01   5.52e-03   1.02e-03s
1600   6.9733e+03   2.01e-02   7.19e+01   5.52e-03   1.15e-03s
1800   7.2828e+03   1.32e-01   5.35e+01   5.52e-03   1.28e-03s
2000   7.0420e+03   5.10e-02   4.84e+01   5.52e-03   1.41e-03s
2200   6.8012e+03   7.67e-02   1.34e+01   5.52e-03   1.53e-03s
2400   6.9411e+03   4.85e-02   5.56e+01   5.52e-03   1.66e-03s
2600   7.2646e+03   1.14e-01   3.67e+01   5.52e-03   1.79e-03s
2800   7.0907e+03   4.50e-02   5.54e+01   5.52e-03   1.92e-03s
3000   6.7481e+03   1.40e-01   1.47e+01   5.52e-03   2.05e-03s
3200   6.9816e+03   3.59e-02   4.20e+01   5.52e-03   2.18e-03s
3400   7.2479e+03   9.89e-02   3.30e+01   5.52e-03   2.30e-03s
3425   7.2426e+03   8.85e-02   1.51e+00   5.52e-03   2.32e-03s
plsh   7.0172e+03   9.01e-14   8.21e-11   --------   2.37e-03s

status:               solved
solution polish:      successful
number of iterations: 3425
optimal objective:    7017.2341
run time:             2.37e-03s
optimal rho estimate: 3.56e-02

Warning: Quadratic matrix (H) is empty or all zero. Your problem is linear. Consider calling LINPROG instead. 
> In quadprog (line 336)
  In dummy (line 20) 

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

max constraint violation = 4.54747e-13
max difference in x = 5.125e-10
obj val difference = 8.93124e-10

Note that the final line of the progress updates is plsh in this output with errors that are all machine precision-ish, and it reports a few lines below that polishing it worked.

Sorry, you’re right. I forgot to delete the opt (which I had been playing with) before posting the results from an empty opt.

And, you’re right that turning on polish does work with this code. Not sure what happened, but it looks like something small must have changed when I attempted to extract an MWE from the environment where I was seeing the failures even with polish on.

I know I had a version of the problem that was giving a different answer, polish was failing, and attempts to tighten bounds only made it terminate at max iterations. I thought I’d confirmed all of that on the exact problem I posted, but apparently not.

Let me have a look and see if I can find what happened and I’ll get back. Sorry about that.

Ok, try this one … a slightly different version of the problem. This one solves in 250 iterations, but polish fails and tightening tolerances also fails.

% problem data
P = sparse(14, 14);
q = [zeros(11,1); 2403.5; 1; 1];
i = [1 1 1 2 2 2 3 3 3 4 4 4 4 5 5 5 6 6 6 6 7 7 7 8 8 8 8 9 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 22 23 24]';
j = [1 4 10 2 8 12 3 6 11 1 4 5 9 4 5 6 3 5 6 7 6 7 8 2 7 8 9 4 8 9 5 6 6 7 7 8 2 8 8 9 4 9 10 13 10 13 10 13 11 14 11 14 1 10 11 12]';
a = [17.3611111 -17.3611111 -1 16 -16 -1 17.0648464 -17.0648464 -1 -17.3611111 39.9953822 -10.8695652 -11.7647059 -10.8695652 16.7519182 -5.88235294 -17.0648464 -5.88235294 32.8678343 -9.92063492 -9.92063492 23.8095238 -13.8888889 -16 -13.8888889 36.100069 -6.21118012 -11.7647059 -6.21118012 17.975886 5.88235294 -5.88235294 9.92063492 -9.92063492 13.8888889 -13.8888889 -16 16 6.21118012 -6.21118012 -11.7647059 11.7647059 2500 -1 3000 -1 3500 -1 1500 -1 2000 -1 1 1 1 1]';
A = sparse(i, j, a, 24, 14);
l = [0 0 0 0 -0.9 0 -1 0 -1.25 -1.5 -0.4 -2.5 -2.5 -2.5 -2.5 -Inf -Inf -Inf -Inf -Inf 0 0.9 0.1 0.1]';
u = [0 0 0 0 -0.9 0 -1 0 -1.25 1.5 0.4 2.5 2.5 2.5 2.5 0 500 1500 0 1000 0 2.5 2.7 3]';

% solve via OSQP
prob = osqp;
prob.setup(P, q, A, l, u, 'polish', 1);
results = prob.solve();
lb_res = A * results.x - l;
ub_res = u - A * results.x;

% solve via quadprog
[x, f, eflag] = quadprog(P, q, [A;-A], [u;-l]);

% compare results
fprintf('max constraint violation = %g\n', abs(min([lb_res; ub_res])))
fprintf('max difference in x = %g\n', norm(x - results.x, Inf));
fprintf('obj val difference = %g\n', norm(f - results.info.obj_val, Inf));

Results in the following for me:

-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 14, constraints m = 24
          nnz(P) + nnz(A) = 56
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, 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: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -6.6980e+03   6.11e+00   7.48e+05   1.00e-01   7.08e-05s
 200   6.7746e+03   1.84e-02   1.86e+03   4.30e-02   1.77e-04s
 250   7.5426e+03   6.31e-01   2.24e-01   6.65e-02   2.09e-04s

status:               solved
solution polish:      unsuccessful
number of iterations: 250
optimal objective:    7542.5719
run time:             2.40e-04s
optimal rho estimate: 3.53e+00

Warning: Quadratic matrix (H) is empty or all zero. Your problem is linear.
Consider calling LINPROG instead. 
> In quadprog (line 336)
  In osqp_lp_fail_2 (line 19) 

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

max constraint violation = 0.0037351
max difference in x = 1280.14
obj val difference = 746.17

I should also mention that it seems that tightening eps_prim_inf and eps_dual_inf seem to have no effect. It still terminates with a solved status and constraint violation > 1e-3. What am I missing?

I was able to get it to solve after playing a bit with the settings, e.g. this works:

prob.setup(P, q, A, l, u, 'polish', 1, 'eps_rel',5e-4, 'rho',0.01,'adaptive_rho_interval',1000,'max_iter',20000,'scaling',1);

It also converged a bit faster when I scaled the linear term in the cost to have unit norm.

However, I would certainly agree that the user should not need to play around so much with the settings to get the solver to converge.

I think the issue in this case is that the “rho adapation” feature just isn’t working very well for this problem, though it is unclear why. A general observation we have had is that the solver tends to work better on QPs than on (some) LPs. If you are willing to post this problem as a github issue we can put it into our bank of awkward problems.

The error terms eps_prim_inf and eps_dual_inf are only used by the solver when trying to identify if the problem is actually infeasible. Your problem is not infeasible, so changing these parameters should not have any effect.

The primal and dual residual terms you are seeing reported are bigger than the setting you gave the solver. This is because the setting you give is the relative tolerance, i.e. primal (or dual) residual scaled by some norm of the problem data. The number reported in the solver progress output is (If I remember right) the absolute error.

Thanks for your help.

I have submitted a GitHub issue with this problem.

And I had just figured out what eps_prim_inf and eps_dual_inf were for by skimming the main OSQP paper right before I got your reply.

Thanks again for your help and for a very nice piece of open-source software. FYI, I plan to include OSQP as one of the solver options in the modeling and optimization layer (MP-Opt-Model) of my open-source electric power system software MATPOWER.