Matlab crashes when calling osqp.codegen() if user sets

settings.scaling = 0.

The crash occurs in the call to osqp_mex('get_workspace')

Example

%small problem
P = sparse([4., 1.; 1., 2.]);
q = [1; 1];
A = sparse([1., 1; 1, 0; 0, 1]);
l = [1.0; 0.0; 0.0];
u = [1.0; 0.7; 0.7];

%1 = ok, 0 = crash
settings.scaling = 0;

% Solve with osqp and export
solver = osqp;
solver.setup(P, q, A, l, u, settings);
solver.solve();

%crashes here
solver.codegen('export', 'project_type','Xcode','parameters', 'matrices', 'mexname', 'emosqp');

Python 3.6.3 (v3.6.3:2c5fed8, Oct 3 2017, 18:11:49) [MSC v.1900 64 bit (AMD64)] on win32
Type "copyright", "credits" or "license()" for more information.

import sys
print(sys.version)
3.6.3 (v3.6.3:2c5fed8, Oct 3 2017, 18:11:49) [MSC v.1900 64 bit (AMD64)]
import osqp
RuntimeError: module compiled against API version 0xc but this version of numpy is 0xb
Traceback (most recent call last):
File "<pyshell#2>", line 1, in
import osqp
File "C:\Users\Joann\AppData\Local\Programs\Python\Python36\lib\site-packages\osqp_init_.py", line 1, in
from osqp.interface import OSQP
File "C:\Users\Joann\AppData\Local\Programs\Python\Python36\lib\site-packages\osqp\interface.py", line 6, in
import osqp._osqp as _osqp # Internal low level module
ImportError: numpy.core.multiarray failed to import
import numpy

I can't install OSQP via pip on linux. This actually makes it impossible to install the new cvxpy 1.0 on linux. I'd really appreciate it if you could upload something to pypi as soon as possible.

I want to minimize an indefinite quadratic function with both equality and inequality constraints that may get violated depending on various factors. So I want to use l1 penalty method that penalizes the violating constraints.

for example,
I have modified an example, to violate the constraints.

import osqp
import scipy.sparse as sparse
import numpy as np

# Define problem data
P = sparse.csc_matrix([[4., 1.], [1., 2.]])
q = np.array([1., 1.])
A = sparse.csc_matrix([[1., 0.], [0., 1.], [1., 0.], [0., 1.]])
l = np.array([0., 0., 0.2, 1.1])
u = np.array([1., 1., 0.2, 1.1])

# Create an OSQP object
prob = osqp.OSQP()

# Setup workspace and change alpha parameter
prob.setup(P, q, A, l, u, alpha=1.0)

# Solve problem
res = prob.solve()
print res.x

Obviously, this is an infeasible problem, so we need to change the objective function to penalize the error. So, I need help to formulate this problem that can be solved using osqp's python interface.

It seems LDL and AMD components of SuiteSparse are packaged within this Apache 2.0 project.

Now without being an expert and while it's possible, that older versions of these sources are not (L)GPL, i do think the current ones are and this looks incompatible to Apache 2.0 to me.

Even more problematic is the fact, that this project does not contain any original license of LDL/AMD at all, despite being mentioned to be looked at in the sources (which are part of this).

ECOS seems to do it right: kept the Readme's; whole project using GPL too (which is annoying for the user of course)

Code:

import numpy
import scipy.sparse
import osqp

Q = numpy.matrix([[1e4, 0], [0, 0]])
P = scipy.sparse.block_diag([scipy.sparse.kron(scipy.sparse.eye(3), Q), scipy.sparse.kron(scipy.sparse.eye(2), 1e-6)])

q = numpy.zeros(8)

Leq = numpy.matrix([[130, -26000,      0.,      0.,      0.,      0.]])
Ueq = Leq

Lineq = numpy.matrix([[-numpy.inf,    -numpy.inf]])
Uineq = numpy.matrix([[numpy.inf,     numpy.inf]])
L = numpy.hstack([Leq, Lineq])
U = numpy.hstack([Ueq, Uineq])

L_cons = L.copy()
L_cons[0, -2] = L_cons[0, -1] = -20000

U_cons = U.copy()
U_cons[0, -2] = U_cons[0, -1] = 20000

AA = numpy.matrix([[1, -6e-5], [4, 1]])
AAeq = scipy.sparse.block_diag([scipy.sparse.kron(scipy.sparse.eye(3), -numpy.eye(2))])
AAeq += scipy.sparse.block_diag([scipy.sparse.kron(scipy.sparse.eye(3, k=-1), AA)])

BB = numpy.matrix([[-2e-9], [8e-5]])
BBeq = scipy.sparse.vstack([numpy.zeros([2, 2]), scipy.sparse.kron(scipy.sparse.eye(2), BB)])

Aeq = scipy.sparse.hstack([AAeq, BBeq])

AAineq = numpy.zeros([2, 6])
BBineq = scipy.sparse.eye(2)

Aineq = scipy.sparse.hstack([AAineq, BBineq])

A = scipy.sparse.vstack([Aeq, Aineq])

prob1 = osqp.OSQP()

prob1.setup(P, q.T, A, L.T, U.T, warm_start=True)

x1 = prob1.solve().x

prob2 = osqp.OSQP()

prob2.setup(P, q.T, A, L_cons.T, U_cons.T, warm_start=True)

x2 = prob2.solve().x

print("\tx1\t\t\t\tx2")
for i in range(len(x1)):
    if abs(x1[i] - x2[i]) < 100:
        print("{:10.2f}\t\t{:10.2f}".format(x1[i], x2[i]))
    else:
        print("{:10.2f}\t\t{:10.2f}\t*".format(x1[i], x2[i]))

print("\n\n* - large differences")

Output:

-----------------------------------------------------------------
problem:  variables n = 8, constraints m = 8
          nnz(P) + nnz(A) = 25
settings: linear system solver = suitesparse ldl,
          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

iter   objective    pri res    dua res    rho        time
   1   0.0000e+00   2.60e+04   1.28e+08   1.00e-01   4.48e-05s
  25   2.5955e+08   8.00e+00   6.05e+02   1.00e-01   6.62e-05s

status:               solved
number of iterations: 25
optimal objective:    259550024.4092
run time:             7.55e-05s
optimal rho estimate: 5.76e-02

-----------------------------------------------------------------
           OSQP v0.2.1  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2017
-----------------------------------------------------------------
problem:  variables n = 8, constraints m = 8
          nnz(P) + nnz(A) = 25
settings: linear system solver = suitesparse ldl,
          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

iter   objective    pri res    dua res    rho        time
   1   0.0000e+00   2.60e+04   1.28e+08   1.00e-01   3.20e-05s
  25   2.5955e+08   8.00e+00   6.05e+02   1.00e-01   5.17e-05s

status:               solved
number of iterations: 25
optimal objective:    259550096.1118
run time:             6.06e-05s
optimal rho estimate: 5.76e-02

    x1                x2
   -129.98           -129.98
  25992.00          25992.00
   -131.55           -131.55
  25471.16          25472.09
   -133.08           -133.08
  24944.75          24945.90
 -11681.23             -4.59    *
  -2661.77             -1.05    *

* - large differences

It would be expected that adding constraints that are wider than the output given by x1 should not affect x2.

Finalize Rust interface from ebarnard/osqp.rs

• Integrate documentation
• Automate deployment

Hi Bart,

I want to solve the following simple LP (using Julia):

min [1;2;3;4]'*[x1;x2;x3;x4]
s.t. x >= 1, x2 >= 5, x1 + x3 >= 4

The following code

using OSQP


c = [1.; 2; 3; 4]
A = sparse([eye(4);0. 1 0 0; 1. 0 1 0])
l = [1.;1;1;1;5;4]
u = Inf*ones(6)
P = spzeros(size(A,2),size(A,2))

m = OSQP.Model()
OSQP.setup!(m; P=P, q=c, A=A, l=l, u=u)
results = OSQP.solve!(m)

gives me a result close to the trivial solution: OSQP.Results([2.99224, 5.00009, 1.00388, 1.0] and optimal objective: 20.0041.

However, if I formulate the problem with equality constraints Ax=b and slack variables instead of inequality constraints l<=Ax<=u:

# reformulated problem with slack variables and equality condition
As = sparse([eye(4) -eye(4) zeros(4,2); 0 1 0 0 zeros(1,4) -1 0; 1 0 1 0 zeros(1,4) 0 -1])
bs = [ones(4);5;4]
cs = [c;zeros(6)]
Ps = spzeros(size(As,2),size(As,2))

m2 = OSQP.Model()
OSQP.setup!(m2; P=Ps, q=cs, A=As, l=bs, u=bs)
results2 = OSQP.solve!(m2)

the solver returns dual infeasible.

Lines 288 and 291 dereference work if PRINTING is defined. Just a few lines later on line 302 the code checks to see if work is NULL. This either suggests work cannot be NULL or there is a null pointer deference. Since work is passed by the caller, this appears to be a possible null pointer dereference.

https://github.com/oxfordcontrol/osqp/blob/88277e9f0a2c7079002bb12356bd5e5363b42042/src/osqp.c#L288
https://github.com/oxfordcontrol/osqp/blob/88277e9f0a2c7079002bb12356bd5e5363b42042/src/osqp.c#L291
https://github.com/oxfordcontrol/osqp/blob/88277e9f0a2c7079002bb12356bd5e5363b42042/src/osqp.c#L302

osqp_demo.zip
Thanks for providing the solver, it's really fast!

I have a question to impose one side of the bound to infinity.

I tried to solve the following problem

min x₁² + 2x₂²
s.t x₁ + 2x₂ = 2
     x₁ ≥  1
     x₂ ≥  2

Note that one side of the bounds are infinity. The code is attached (I just modified a few entries in examples/osqp_demo.c). The solver reaches iteratation limit, and keeps increasing the weight on the penalty term.

If I change the last two constraints to be upper-bounded by finite numbers, such as

2 ≥ x₁ ≥  1
3 ≥ x₂ ≥  2

then the solver detects the infeasibility immediately.

How should I impose constraint with one side of the bounds being infinity?

Hello and congratulations on this very nice software package! I followed the instructions on how to install the Julia interface to the solver, but it looks like OSQP.jl is not yet a registered package.
I suppose the correct procedure should be that of issuing

Pkg.clone("https://github.com/oxfordcontrol/OSQP.jl.git")

rather than using Pkg.add.

I tried to solve a problem with polish on. The polish fails, and valgrind tells me there is a memory leak.

The test code is attached, it is a modification of osqp_demo.c.

To run the test, I set DFLOAT to ON in cmake, and ran

$ cd OSQP_DIRECTORY/build
$ make
$ valgrind --tool=memcheck --leak-check=full out/osqp_demo

The valgrind error message is

-----------------------------------------------------------------
           OSQP v0.2.1  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2017
-----------------------------------------------------------------
problem:  variables n = 3, constraints m = 7
          nnz(P) + nnz(A) = 20
settings: linear system solver = suitesparse ldl,
          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

iter   objective    pri res    dua res    rho        time
   1   0.0000e+00   1.00e+00   1.00e+02   1.00e-01   7.30e-02s
  25   0.0000e+00   3.94e-04   1.21e-04   1.00e-01   8.04e-02s

status:               solved
solution polish:      unsuccessful
number of iterations: 25
optimal objective:    0.0000
run time:             8.53e-02s
optimal rho estimate: 1.98e-03

==30098==
==30098== HEAP SUMMARY:
==30098==     in use at exit: 356 bytes in 4 blocks
==30098==   total heap usage: 136 allocs, 132 frees, 11,460 bytes allocated
==30098==
==30098== 356 (56 direct, 300 indirect) bytes in 1 blocks are definitely lost in loss record 4 of 4
==30098==    at 0x4C2FB55: calloc (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==30098==    by 0x403669: csc_calloc (cs.c:13)
==30098==    by 0x403743: csc_spalloc (cs.c:37)
==30098==    by 0x403F2D: csc_symperm (cs.c:151)
==30098==    by 0x40F6AC: permute_KKT (suitesparse_ldl.c:125)
==30098==    by 0x40FA04: init_linsys_solver_suitesparse_ldl (suitesparse_ldl.c:209)
==30098==    by 0x40F0EC: init_linsys_solver (lin_sys.c:56)
==30098==    by 0x409C4A: polish (polish.c:212)
==30098==    by 0x407164: osqp_solve (osqp.c:464)
==30098==    by 0x4035B5: main (osqp_demo.c:49)
==30098==
==30098== LEAK SUMMARY:
==30098==    definitely lost: 56 bytes in 1 blocks
==30098==    indirectly lost: 300 bytes in 3 blocks
==30098==      possibly lost: 0 bytes in 0 blocks
==30098==    still reachable: 0 bytes in 0 blocks
==30098==         suppressed: 0 bytes in 0 blocks
==30098==
==30098== For counts of detected and suppressed errors, rerun with: -v
==30098== ERROR SUMMARY: 1 errors from 1 contexts (suppressed: 0 from 0)

I think KKT_temp allocated in suitesparse_ldl.c:125 is not de-allocated. It is supposed to be de-allocated in https://github.com/oxfordcontrol/osqp/blob/master/lin_sys/direct/suitesparse/suitesparse_ldl.c#L251, but it never enters that if (polish) branch.

I do not see this memory leaking problem when polish succeeds. For example, if I set DFLOAT to OFF, the polish succeeds, and no memory leak.
osqp_demo.zip

It is useful in soft real-time mpc to be able to have large optimisations run in the free time between mpc iterations. For this to work it must be possible to interrupt the solver at arbitrary times.

This is almost possible using check_termination with a value of 1 and the ctrl-c handler. However, this does not work on windows as signal handlers are per-process so this breaks with multiple osqp solvers running (and does not restore the current signal handler).

I'd suggest adding an interrupt callback that executes at the same time the ctrl-c handler and can be enabled irrespective of the value of the CTRLC build option using:

c_int osqp_set_interrupt_callback(OSQPWorkspace* work, c_int (*callback)(void* opaque), void* opaque);

The callback is disabled if a value of NULL is passed to callback.

I'm happy to implement this. Just want feedback.

Python interface gives segfault when trying to solve problem CONT-200.

Seems to crash while polishing for problems with no active constraints at optimality.

Example:
A = sparse(eye(2));
P = sparse(eye(2));
lA = [-2;-2];
q = [-3;-4];
uA = [10;10];

Probably polishing should just terminate if mred = 0 in polish(). Should also test what happens when polishing runs on a solution that is very close to the unconstrained optimizer.

I have seen in your C/C++ api that it is possible to update both the hessian matrix Pand the linear constraint matrix A using respectively the osqp_update_P() and the osqp_update_A() methods; however these functions assume that the sparsity structure of the matrices does not change.

If this assumption is no more satisfied is it possible to remove the old matrices and add the new ones?
I have seen in your repository that the linsys_solver struct depens on A and P matrices and more specifically the KKT matrix depends on the number of non zero entries in the A and P matrices (so the linsys_solver->update_matrices() method cannot be used for this purpose). Thus my idea is to deallocate the linsys_solver struct and instantiate it again with the new matrices.

...

// Free linear system solver structure
if (work->linsys_solver){
   if(work->linsys_solver->free){
      work->linsys_solver->free(work->linsys_solver);
    }
}

// Free old hessian and linear constraint matrices
csc_spfree(work->data->P);
csc_spfree(work->data->A);

// Allocate new matrices
work->data->P = csc_to_triu(newP); 
work->data->A = copy_csc_mat(newA);

// Instantiate the solver 
work->linsys_solver = init_linsys_solver(work->data->P, work->data->A,
		                         work->settings->sigma, work->rho_vec,
			                 work->settings->linsys_solver, 0);
...

What do you think about this approach? Is there a better way to solve this problem?
Thank you in advance for any help you can provide.

I'm not too familiar with the settings and how they affect the optimization problem, however when I modify osqp_demo_direct.c with the following values, the solver doesn't converge:

    // Load problem data
    c_float P_x[1] = {2.00004};
    c_int P_nnz = 1;
    c_int P_i[1] = {0};
    c_int P_p[5] = {0, 1, 1, 1, 1};
    c_float q[4] = {-2.00035, 0, 1, 1};
    c_float A_x[8] = {-200, 1, 10, 1, 1, 1, -1, 1};
    c_int A_nnz = 8;
    c_int A_i[8] = {0, 1, 0, 2, 0, 3, 0, 4};
    c_int A_p[5] = {0, 2, 4, 6, 8};
    c_float l[5] = {-1000, 9.9, 0.9, 0, 0};
    c_float u[5] = {-1000, 10.1, 1.1, 1e+30, 1e+30};
    c_int n = 4;
    c_int m = 1;

Currently the Python API requires P and A to be in CSC format.

Problem is: on small instances, it can take much more time to convert to CSC than to solve the QP itself. In this small example where P and A are 3x3, CSC conversions takes ~200 µs on my machine, versus ~20 µs taken to solve the QP itself.

In [17]: P.shape
Out[17]: (3, 3)

In [18]: %timeit csc_matrix(P)
10000 loops, best of 3: 116 µs per loop

Could the API allow for dense numpy.array arguments for P and A?

Code:

import osqp

prob = osqp.OSQP()
prob.update_settings(verbose=False)

Output:

Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)

Did not expect error.

Ran on OSQP 0.2.1 with Python 3.6.4 installed via Anaconda on Ubuntu 16.04.
Replicated for OSQP 0.2.1 with Python 3.6.4 installed via Anaconda on Mac OS X High Sierra.

OSQP currently enforces that eps_rel be positive:
https://github.com/oxfordcontrol/osqp/blob/515ceba7c72c47ea73f01b136cb729d62ae27dc0/src/auxil.c#L887-L890

Is there a reason to not allow eps_abs positive and eps_rel zero? Many other QP solvers allow you to enforce an absolute tolerance on constraints, so setting eps_rel to zero would help when comparing OSQP to these solvers. Of course, it's possible to set it to an extremely tiny value, but you shouldn't need to do this.

When an optimization problem is initialized all the data saved in the OSQPData struct are copied inside OSQPWorkspace by the way the number of non-zero elements of the hessian matrix seems to be non-coherent.

As regarding the Hessian matrix P only its upper-diagonal part is saved by using
https://github.com/oxfordcontrol/osqp/blob/71aaadb3daad3ae133400d1d0b73f1644cd21ae9/src/osqp.c#L114
Inside this method the number of non-zero elements is estimated
https://github.com/oxfordcontrol/osqp/blob/71aaadb3daad3ae133400d1d0b73f1644cd21ae9/src/cs.c#L259-L266
and then saved
https://github.com/oxfordcontrol/osqp/blob/71aaadb3daad3ae133400d1d0b73f1644cd21ae9/src/cs.c#L309-L310

However the numbers of non-zero elements of data->P and work->data->P seems to be different.
I wrote a simple test in order to better explain this behaviour osqp-test.zip.
In the following some lines are shown

// Structures
OSQPWorkspace * work;  // Workspace
OSQPData * data;  // OSQPData

// Populate data
data = (OSQPData *)c_malloc(sizeof(OSQPData));
data->n = n;
data->m = m;
data->P = csc_matrix(data->n, data->n, P_nnz, P_x, P_i, P_p);
data->q = q;
data->A = csc_matrix(data->m, data->n, A_nnz, A_x, A_i, A_p);
data->l = l;
data->u = u;

// Define Solver settings as default
osqp_set_default_settings(settings);

settings->alpha = 1.0; // Change alpha parameter

// Setup workspace
work = osqp_setup(data, settings);
printf("#nonzero elements (original) %lld\n", data->P->nzmax);
printf("#nonzero elements (inside workspace) %lld\n", work->data->P->nzmax);

The output is

-----------------------------------------------------------------
           OSQP v0.3.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2017
-----------------------------------------------------------------
problem:  variables n = 2, constraints m = 3
          nnz(P) + nnz(A) = 8
settings: linear system solver = suitesparse ldl,
          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.00, max_iter = 4000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: off

#nonzero elements (original) 2
#nonzero elements (inside workspace) 4

Is this behaviour desired? Is it a bug? If it's a bug I think it can be easily solved by changing this line
https://github.com/oxfordcontrol/osqp/blob/71aaadb3daad3ae133400d1d0b73f1644cd21ae9/src/cs.c#L309-L310
into

M_triu->nzmax = z_M; 

Solving the unconstrained problem of the form

minimize   1/2 x' P x + q'x

does not converge even though primal and dual residuals are both zero during the iterations.

OSQP has some local modifications to LDL. I think the diff is limited to d282e7a. Given that the modifications are trivial (rescaling a vector), could they be separated from the local copy of LDL? I'm in a situation where I'd like to link OSQP against vanilla SuiteSparse.

Polishing procedure fails when the eps_abs and eps_rel tolerances are below the reduced KKT regularization term; i.e. the polished solution is less accurate than the one obtained by the ADMM algorithm.

We are using the following perturbed reduced KKT matrix

KKTred + [delta*eye 0; 0 -delta*eye]

to enforce quasi-definiteness.

• I think we should apply iterative refinement for a few steps like in ECOS (Sec IV-A.4) so that the polished solution accuracy increases.

• There should be anyway a presolve-check to see whether it makes sense to polish if the defined tolerances are already low. If the user sets very low tolerances, it could be useful to stop the ADMM algorithm at certain tolerances (different than the user ones) and polish to obtain a more accurate solution.

An error occurs when building the source files on macOS (10.12.5).

[ 68%] Generating unittests data files using Python
Traceback (most recent call last):
File "generate_tests_data.py", line 3, in
import basic_qp.generate_problem
File "/Users/user/Desktop/OSQP_TEST/osqp/tests/basic_qp/generate_problem.py", line 3, in
import utils.codegen_utils as cu
File "/Users/user/Desktop/OSQP_TEST/osqp/tests/utils/codegen_utils.py", line 3, in
from builtins import dict
ImportError: No module named builtins
make[2]: *** [../tests/basic_qp/data.h] Error 1
make[1]: *** [CMakeFiles/osqp_tester_direct.dir/all] Error 2
make: *** [all] Error 2

The macOS (10.12.5) has a built-in Python 2.7 library, and CMAKE will use the built-in Python library if not specified. The module builtins is a Python3 module, and it has a different name __builtin__ in Python2.

My solution is to add a Python version requirement in CMakeLists.txt, and then CMAKE will find the correct Python3 library I installed.

It should be also noticed that the function random in scipy.sparse (used in /osqp/tests/utils/codegen_utils.py) was added in scipy after v0.17.0, and there should also be a version requirement for scipy.

This line:
https://github.com/oxfordcontrol/osqp/blob/2aeaebb3a3f8f718e16373caccb0c5a8e091a004/src/osqp.c#L343
will attempt to convert NaN to a c_int in the case that work->settings->check_termination equals zero, which is a valid value according to:
https://github.com/oxfordcontrol/osqp/blob/2aeaebb3a3f8f718e16373caccb0c5a8e091a004/include/types.h#L167

This conversion from NaN to integer is required to trigger an invalid operation exception according to the IEEE 754-2008 standard. On most systems this is silently ignored, but on others, under certain compiler settings, this causes the program to crash.

I was able to work around the issue by wrapping the offending line with if (work->settings->check_termination > 0) {}, but I don't understand the context of the code to know if this is correct.

If the number of scaling iteration is set to 0 a segmentation fault occurs when the hessian P or/and the linear constraints matrix A are updated.
You can verify the issue adding the following line here

settings->scaling = 0;

I think that the problem is caused to the calling of the scale_data() and unscale_data() methods in the osqp_update_*() functions. In this case it can be easily solved adding an if statement

if(work->settings->scaling){
    scale_data(work);
}

This is a feature request.

Would be great to use MUMPS (http://mumps.enseeiht.fr) as the direct solver in OSQP.

When OSQP is compiled with DLONG turned on and I link to it from my project, the project crashes. However when I explicitly add "add_definitions(-DDLONG)" to my cmake file, everything works fine. This is a problem because once the user has compiled OSQP there's no way to tell what options they've compiled with.

I tried to write some function that will link to both SCS and OSQP. Unfortunately SCS and OSQP share some functions with the same signature, that SCS calls OSQP functions internally, and resulted in a segmentation fault. Here is the valgrind output.

==25464== Use of uninitialised value of size 8
==25464==    at 0x40339CA: ldl_l_perm (ldl.c:424)
==25464==    by 0x4034A70: LDLSolve (suitesparse_ldl.c:284)
==25464==    by 0x4050D77: scs_solve_lin_sys (private.c:224)
==25464==    by 0x405522C: update_work (scs.c:820)
==25464==    by 0x4054BAF: scs_solve (scs.c:853)
==25464==    by 0x40576E1: scs (scs.c:979)

Note that in private.c:224, it calls the LDLSolve function defined in OSQP https://github.com/oxfordcontrol/osqp/blob/master/lin_sys/direct/suitesparse/suitesparse_ldl.c#L280~L291. Do you think you can change the CMakeLists.txt file, and use the trick in https://cmake.org/cmake/help/v3.5/module/GenerateExportHeader.html to hide this function?

I also submitted an issue to SCS in cvxgrp/scs#94

I get the following in from the iPython console:
An error ocurred while starting the kernel In file included from ../include/types.h:9:0, from ../include/osqp.h:9, from emosqpmodule.c:7: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/ldl.h:10, from osqp/ldl.c:157: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/osqp.h:9, from osqp/auxil.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/scaling.h:9, from osqp/scaling.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/proj.h:8, from osqp/proj.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/osqp.h:9, from osqp/osqp.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/ldl.h:10, from osqp/suitesparse_ldl.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/lin_alg.h:9, from osqp/lin_alg.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/util.h:8, from osqp/util.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^ In file included from ../include/types.h:9:0, from ../include/kkt.h:8, from osqp/kkt.c:1: ../include/constants.h:34:20: warning: ‘LINSYS_SOLVER_NAME’ defined but not used [‑Wunused‑variable] static const char *LINSYS_SOLVER_NAME[] = { ^

I suggest that we indicate norm types, in equilibration and elsewhere, via a custom type. For example, we should define something like

typedef enum {ONE_NORM, TWO_NORM, INF_NORM} norm_t;

and require use of this type when performing equilibration. The current (1,2,-1) integer flags are not that nice, and possible error prone.

Running the 'run_osqp_tests.m' scripts crashes in R2016a.

The crash happens during the codegen tests, since the function osqp.codegen expects to get a 'folder' field from the matlab 'dir' command (see osqp.m line 425). This behaviour of the matlab 'dir' command seems to have been added from R2016b.

The matlab interface should implement osqp.default_settings as a static method.

Feasibility problem of the form

find            0
subject to  A x = b

converges to the right solution, but the resulting status is "Infeasible".

New custom LDL solver for OSQP.

osqp failed to find a feasible solution when solving a quadratic programming problem given below. Even if I provide an initial feasible solution using prob.warm_start(x=x0), it will still fail to find a solution on the first try. However, the second try will succeed if I use prob.warm_start(x=x0) after the first try.

It seems to me that I have to set the value of the y correctly before calling solve(). However, it is not clear to me how to do this correctly.

python script to reproduce the error:

import pandas as pd
import scipy as sp
import numpy as np
import osqp
from cvxopt import solvers, matrix
P, q, A, l, u, x0 = pd.read_msgpack('args.msgpack')

prob = osqp.OSQP()
prob.setup(P=sp.sparse.csc_matrix(P), q=q, A=sp.sparse.csc_matrix(A), l=l, u=u)
prob.warm_start(x=x0)
res1 = prob.solve()
prob.warm_start(x=x0)
res2 = prob.solve()

res3 = solvers.qp(matrix(P), matrix(q), matrix(A), matrix(u))

print('')
print('Since A x < u, min(u - Ax) > 0.')
print('x0 : min(u - A x) = {}.'.format((u - A @ x0).min()))
print('osqp_warm1 : min(u - A x) = {}.'.format((u - A @ res1.x).min()))
print('osqp_warm2 : min(u - A x) = {}.'.format((u - A @ res2.x).min()))
print('cvxopt : min(u - A x) = {}.'.format((u - A @ np.array(res3['x'])[:, 0]).min()))

output:

-----------------------------------------------------------------
           OSQP v0.2.1  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2017
-----------------------------------------------------------------
problem:  variables n = 450, constraints m = 755
          nnz(P) + nnz(A) = 6008
settings: linear system solver = suitesparse ldl,
          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

iter   objective    pri res    dua res    rho        time
   1  -1.0511e+00   7.30e+00   4.00e-04   1.00e-01   9.48e-04s
 200  -2.6917e+30   4.36e+17   2.10e+15   2.78e+00   5.04e-03s
 400  -8.4202e+88   5.92e+47   2.50e+39   1.00e+06   9.69e-03s
 600  -8.4013e+88   3.27e+47   4.52e+39   1.00e+06   1.37e-02s
 800  -8.3824e+88   5.92e+47   2.87e+39   1.00e+06   1.76e-02s
1000  -8.3636e+88   4.37e+47   3.31e+39   1.00e+06   2.16e-02s
1200  -8.3448e+88   4.73e+47   3.20e+39   1.00e+06   2.55e-02s
1400  -8.3261e+88   6.49e+47   2.59e+39   1.00e+06   2.94e-02s
1600  -8.3074e+88   4.52e+47   3.07e+39   1.00e+06   3.34e-02s
1800  -8.2887e+88   6.52e+47   3.20e+39   1.00e+06   3.74e-02s
2000  -8.2701e+88   9.19e+47   3.44e+39   1.00e+06   4.17e-02s
2200  -8.2516e+88   8.45e+47   2.54e+39   1.00e+06   4.59e-02s
2400  -8.2331e+88   9.22e+47   3.05e+39   1.00e+06   5.01e-02s
2600  -8.2146e+88   2.86e+47   2.29e+39   1.00e+06   5.42e-02s
2800  -8.1962e+88   6.54e+47   2.44e+39   1.00e+06   5.84e-02s
3000  -8.1778e+88   6.83e+47   3.48e+39   1.00e+06   6.26e-02s
3200  -8.1594e+88   4.75e+47   4.00e+39   1.00e+06   6.69e-02s
3400  -8.1411e+88   3.26e+47   3.01e+39   1.00e+06   7.11e-02s
3600  -8.1229e+88   6.50e+47   3.80e+39   1.00e+06   7.53e-02s
3800  -8.1047e+88   2.97e+47   2.28e+39   1.00e+06   7.95e-02s
4000  -8.0865e+88   6.48e+47   2.81e+39   1.00e+06   8.37e-02s

status:               maximum iterations reached
number of iterations: 4000
run time:             8.37e-02s
optimal rho estimate: 1.00e+06

iter   objective    pri res    dua res    rho        time
   1  -1.2255e-07   1.33e-01   1.33e+02   1.00e+06   5.58e-05s
  25   1.5934e-03   1.06e-06   1.00e-03   1.00e+06   5.44e-04s

status:               solved
number of iterations: 25
optimal objective:    0.0016
run time:             5.84e-04s
optimal rho estimate: 2.59e+03

     pcost       dcost       gap    pres   dres
 0: -1.6528e-04 -2.8047e+01  9e+02  1e+01  3e+01
 1:  1.5753e-01 -3.4792e+01  5e+01  3e-01  8e-01
 2:  1.5876e-01 -3.1946e+00  3e+00  1e-15  1e-14
 3:  1.1445e-02 -2.1430e-01  2e-01  9e-16  1e-15
 4:  3.4743e-03 -1.1573e-02  2e-02  2e-16  9e-17
 5:  2.2934e-03 -4.6914e-03  7e-03  1e-16  4e-17
 6:  1.4080e-03 -2.8373e-03  4e-03  1e-16  3e-17
 7:  6.2010e-04 -1.8405e-03  2e-03  1e-16  3e-17
 8:  3.2233e-04 -1.3853e-03  2e-03  1e-16  2e-17
 9:  1.3188e-06 -9.7237e-04  1e-03  1e-16  2e-17
10: -2.4953e-04 -6.5637e-04  4e-04  2e-16  2e-17
11: -2.7798e-04 -6.0490e-04  3e-04  1e-16  3e-17
12: -3.8672e-04 -5.0641e-04  1e-04  1e-16  3e-17
13: -4.4245e-04 -4.6236e-04  2e-05  1e-16  1e-16
14: -4.4564e-04 -4.5654e-04  1e-05  1e-16  4e-16
15: -4.4886e-04 -4.5359e-04  5e-06  1e-16  4e-16
16: -4.5101e-04 -4.5137e-04  4e-07  2e-16  7e-16
17: -4.5114e-04 -4.5115e-04  4e-09  2e-16  1e-15
Optimal solution found.

Since A x < u, min(u - Ax) > 0.
x0 : min(u - A x) = 0.01.
osqp_warm1 : min(u - A x) = -3.905538537264496e+47.
osqp_warm2 : min(u - A x) = -9.560234414311724e-08.
cvxopt : min(u - A x) = 4.72754359270855e-09.

Codes and data are attached in the following zipfile:
osqp_fail.zip

I was a bit surprised to discover that osqp_warm_start_x clears the dual warm start and osqp_warm_start_y clears the primal warm start:
https://github.com/oxfordcontrol/osqp/blob/370dd33622aa1fec1baae3881d7a5a56bf0ebb19/src/osqp.c#L899-L900

https://github.com/oxfordcontrol/osqp/blob/370dd33622aa1fec1baae3881d7a5a56bf0ebb19/src/osqp.c#L918-L920

This isn't documented, and it's natural to assume that these functions would just act like setters for the primal and dual warm starts.

I'm following the instructions:

!git clone --recurse-submodules https://github.com/oxfordcontrol/osqp-matlab
cd osqp-matlab
make_osqp

But I get the following output error:

Error using mex
'/Users/ajx/Dropbox/l1/osqp-matlab/osqp_mex.mexmaci64' is compiled
with incompatible options '-R2017b' and -R2018a'. For more
information, see MEX file compiled with incompatible options.
Error in make_osqp (line 178)
    eval(cmd); 

I'm on Mac OS X 10.12 with MATLAB R2018a

Develop new linear system solver on CUDA GPU.

OSQP performance is great! One thing that would be nice to have: Can it support non-linear user-defined functions directly? Like Julia do?

Examples, constrains like: o[1] == p[1] * (ti[1]^p[4]) + p[2] * ti[1] + p[3]

So I've been trying to use the Python script below to generate the OSQP code for my embedded system, but I've found that it has not been generating arrays of the correct size. More specifically, when I go into the generated workspace.h file, instead of having a P matrix of size 16 (4x4), it is only size 10. I am worried that these incorrect array sizes are being carried through the rest of workspace.h and possibly other parts of the generated code and are causing the solver to fail as a result of that.

# get an OSQP object
m = osqp.OSQP()

# linear cost term
q = np.zeros((1, 4))[0] + 1
# lower wheel acceleration constraint
l = np.zeros((1, 4))[0] + 1
# upper wheel acceleration constraint
u = np.zeros((1, 4))[0] + 1
# quadratic term
P = sparse.csc_matrix(np.zeros((4, 4)) + 1)
# A matrix for l <= Ax <= u
A = sparse.csc_matrix(np.zeros((4, 4)) + 1)

# setup the solver
m.setup(P=P, A=A, q=q, l=l, u=u)
# generate the C code in the folder that this file is run from.
m.codegen(os.getcwd(), parameters="matrices")

Hi,

The problem that I have faced is that when I put a quadratic programming optimization in a loop (say for 1000 runs), it does not come back with results for most of the cases. I traced the problem and kept the record of the failures. Surprisingly, when I run the optimization for those failed ones individually, it does not fail, but it fails when in a loop. Also, the failures are not reproducible, meaning that the loop each time come back with different failed cases!

Hi there, I've been receiving the error

"Could NOT find Matlab (missing: Matlab_MEX_LIBRARY) (found version "9.0")
CMake Error at CMakeLists.txt:220 (message):
You need Matlab libraries to build the Matlab interface"

when building the Matlab interface (using make_osqp command in Matlab). It seems like CMake cannot find Matlab. Any idea to fix this? Thanks.

Attempting to warm start OSQP when scaling is disabled causes a segfault. Example via the matlab interface:

%make a solver with default settings
%and disasble scaling
solver = osqp;
osqpOptions = solver.default_settings();
osqpOptions.scaling = 0;   %uncomment to crash

%load any problem
file = 'problems/maros/marosHS76.mat';
problem = readProblem(file);

%assign data and configure the initial iterates
solver.setup(problem.P,problem.q,problem.A,problem.l,problem.u,osqpOptions);

%configure for warm starting
x = zeros(size(problem.A,2),1);
y = zeros(size(problem.A,1),1);
solver.warm_start('x',x,'y',y);  %crash!

%solve
solver.solve();

Could you change the line below so it matches ECOS?
https://github.com/oxfordcontrol/osqp/blob/df90dafe0c4f1a1ca17f9d04865afbf103d77e1e/include/glob_opts.h#L149

# ECOS version
#elif defined PYTHON
#include <Python.h>
#define PRINTTEXT PySys_WriteStdout

This is important for my unit tests.

On the following example:

import numpy
import scipy.sparse
from osqp import OSQP

if __name__ == "__main__":
    n = 4
    P = scipy.sparse.lil_matrix(scipy.sparse.eye(n))
    for i in xrange(1, n - 1):
        P[i, i + 1] = -1
        P[i, i - 1] = 1
    P = scipy.sparse.csc_matrix(P)
    q = -numpy.ones((n,))
    A = scipy.sparse.csc_matrix(-scipy.sparse.eye(n))
    u = -2 * numpy.ones((n,))
    l = -numpy.inf * numpy.ones(n)
    osqp = OSQP()
    osqp.setup(P=P, q=q, A=A, l=l, u=u)
    res = osqp.solve()
    print "x =", res.x

The expected solution is x = [2 2 2 2], but OSQP returns x = [2 nan nan nan](with a successful return status) on my machine. Here is the full log:

------------------------------------------------------------
        OSQP v0.1.1  -  Operator Splitting QP Solver
           (c) Bartolomeo Stellato,  Goran Banjac
     University of Oxford  -  Stanford University 2017
------------------------------------------------------------
Problem:  variables n = 4, constraints m = 4
Settings: 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 
          sigma = 1.0e-06, alpha = 1.6e+00, 
          max_iter = 2500
          early_terminate: on (interval 25)
          scaling: on, scaled_termination: off
          warm start: on, polish: on

Iter    Obj  Val     Pri  Res     Dua  Res       Time 
   1  -1.2197e+00   8.2785e+00   4.0001e-01      0.00s
 100  -5.5906e+36   3.6738e+18   8.8694e+11      0.00s
 200  -3.1426e+72   2.7544e+36   6.6498e+29      0.00s
 300 -1.7665e+108   2.0651e+54   4.9857e+47      0.00s
 400 -9.9301e+143   1.5483e+72   3.7380e+65      0.00s
 500 -5.5819e+179   1.1609e+90   2.8026e+83      0.00s
 600 -3.1377e+215  8.7035e+107  2.1012e+101      0.00s
 700 -1.7638e+251  6.5255e+125  1.5754e+119      0.00s
 800 -9.9147e+286  4.8925e+143  1.1811e+137      0.00s
 900         -nan  3.6681e+161  8.8556e+154      0.00s
1000         -nan  2.7502e+179  6.6395e+172      0.00s
1100         -nan  2.0619e+197  4.9779e+190      0.00s
1200         -nan  1.5459e+215  3.7322e+208      0.00s
1300         -nan  1.1591e+233  2.7982e+226      0.00s
1400         -nan  8.6900e+250  2.0980e+244      0.00s
1500         -nan  6.5153e+268  1.5729e+262      0.00s
1600         -nan  4.8849e+286  1.1793e+280      0.00s
1700         -nan  3.6624e+304  8.8419e+297      0.00s
1725         -nan   1.5543e-15   4.4409e-16      0.00s

Status: Solved
Solution polish: Unsuccessful
Number of iterations: 1725
Optimal objective: -nan
Run time: 0.264ms

set_default_settings is exported from osqp.h but is missing the osqp_ prefix. This name is pretty generic, is there a reason not to prefix it?
Ideally everything defined by including osqp.h should have a prefix, but set_default_settings is maybe the most prominent candidate because it's part of the public API and is used everywhere.

The osqp.m interface does not allow true/false values as settings. For example, it insists on

>> options.polish = 0;

instead of allowing also

>> options.polish = false;

The adaptive_rho_fraction field is only defined with PROFILING:
https://github.com/oxfordcontrol/osqp/blob/2aeaebb3a3f8f718e16373caccb0c5a8e091a004/include/types.h#L146
but it's accessed here:
https://github.com/oxfordcontrol/osqp/blob/2aeaebb3a3f8f718e16373caccb0c5a8e091a004/src/auxil.c#L845
regardless of PROFILING.
This is with /* #undef EMBEDDED */ in glob_opts.h.