Add QP solver through OSQP

[cdiener]

This is a first working version of a QP solver via OSQP. OSQP expects the problem to be sparse and static so I added a thin compatibility layer that manages the OSQP problem via dictionaries. All operations on that layer are at worst O(nv + nc) but it makes problem setup about 10-20% of the entire solve time for most models. Could maybe be optimized more in the future but should be sufficient for now. Still need to add tests but I added some more complicated models to the test data for now.

It now supports all features such as reduced costs etc. This version is also faster than the previous PR since it gets primals and duals via the C interface instead of evaluating expressions. We used it in the recent MICOM course and it worked well.

Demo:

Python 3.9.1 (default, Dec  8 2020, 07:51:42) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.19.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: from optlang import osqp_interface as o

In [2]: import json

In [3]: qp = json.load(open("src/optlang/tests/data/QPLIB_8785.json"))

In [4]: mod = o.Model.from_json(qp)

In [6]: len(mod.variables), len(mod.constraints)
Out[6]: (10399, 11362)

In [7]: mod.configuration.verbosity = 1

In [8]: mod.configuration.tolerances.optimality = 1e-6

In [9]: mod.optimize()
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 10399, constraints m = 21761
          nnz(P) + nnz(A) = 73422
settings: linear system solver = qdldl,
          eps_abs = 1.0e-06, eps_rel = 1.0e-06,
          eps_prim_inf = 1.0e-06, eps_dual_inf = 1.0e-06,
          rho = 1.00e+00 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 100000
          check_termination: on (interval 25),
          scaling: off, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1   0.0000e+00   9.90e+01   9.90e+01   1.00e+00   5.97e-01s
 200   7.8676e+03   3.39e-02   2.37e-05   9.64e-04   3.62e+00s
 300   7.8675e+03   1.17e-04   2.22e-07   6.19e-03   4.87e+00s
plsh   7.8675e+03   8.53e-14   6.68e-10   --------   5.60e+00s

status:               solved
solution polish:      successful
number of iterations: 300
optimal objective:    7867.4911
run time:             5.60e+00s
optimal rho estimate: 5.22e-03

Out[9]: 'optimal'

In [10]: mod.objective.value  # reference solution: 7867.4911490000004051
Out[10]: 7867.491148760153

One should note that this is a direct solver so it performs pretty bad for LPs or at very low tolerances. I added some optimizations for LP so it works ok-ish but worse than GLPK. However it is great for QPs so I think having a specialized solver for that is still helpful.

TODO

• [x] add tests
• [ ] add documentation
• [x] update deps?

[cdiener]

@Midnighter @KristianJensen this is done for now. There may be small fixes we can make in the future but it's pretty complete.

[KristianJensen]

@cdiener I will try to have a look at it as soon as possible

[cdiener]

Failing now due to missing packages for swiglpk 5.0.0 for Mac and some security issue with pyyaml which is not a dependency so no idea where that comes from.

[cdiener]

Okay nice, at that point the only thing missing is to remove Python 3.5 from the CI. I'll send a separate PR for that.

[cdiener]

Ups, that was a bit too early. This needed additional fixes to make the CI work. I'll send another PR.