[WIP]: Implementation of OSQP algorithm in Proxsuite

.github/workflows/conda/environment.yml

@@ -18,3 +18,4 @@ dependencies:
   - libmatio
   - numpy
   - scipy
+  - osqp

CHANGELOG.md

@@ -8,6 +8,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 ### Added
 - Recursive stub generation for Python bindings ([#419](https://github.com/Simple-Robotics/proxsuite/pull/419))
+- Add OSQP solver and refactor code ([#415](https://github.com/Simple-Robotics/proxsuite/pull/415))
 
 ### Changed
 - Change the default branch to `devel` ([#395](https://github.com/Simple-Robotics/proxsuite/pull/395))

README.md

@@ -166,6 +166,13 @@ where $x \in \mathbb{R}^n$ is the optimization variable. The objective function
 ### Citing **QPLayer**
 
 If you are using **QPLayer** for your work, we encourage you to [cite the related paper](https://inria.hal.science/hal-04133055v2/).
+## **WIP: OSQP**
+
+The **OSQP** algorithm is a numerical optimization approach for solving quadratic programming problems with the same form as problems treated by **ProxQP**. It is based on the Alternating Direction Method of Multipliers.
+
+### Citing **OSQP**
+
+**OSQP** was developped by B. Stellato, G. Banjac, P. Goulart, A. Bemporad and S. Boyd. Information about the algorithm and the API of **OSQP** are available in the [related paper](https://web.stanford.edu/~boyd/papers/pdf/osqp.pdf) and [related website](https://osqp.org/).
 ## Installation procedure
 
 Please follow the installation procedure [here](https://github.com/Simple-Robotics/proxsuite/blob/devel/doc/5-installation.md).

benchmark/timings-box-constraints.cpp

@@ -3,13 +3,13 @@
 //
 #include <iostream>
 #include <proxsuite/proxqp/dense/dense.hpp>
-#include <proxsuite/proxqp/utils/random_qp_problems.hpp>
+#include <proxsuite/common/utils/random_qp_problems.hpp>
 
 using T = double;
 using I = long long;
 
 using namespace proxsuite;
-using namespace proxsuite::proxqp;
+using namespace proxsuite::common;
 
 int
 main(int /*argc*/, const char** /*argv*/)
@@ -20,32 +20,32 @@ main(int /*argc*/, const char** /*argv*/)
   T sparsity_factor = 0.75;
   T eps_abs = T(1e-9);
   T elapsed_time = 0.0;
-  proxqp::utils::rand::set_seed(1);
+  common::utils::rand::set_seed(1);
   std::cout << "Dense QP" << std::endl;
-  for (proxqp::isize dim = 100; dim <= 1000; dim = dim + 100) {
+  for (isize dim = 100; dim <= 1000; dim = dim + 100) {
 
-    proxqp::isize n_eq(dim / 2);
-    proxqp::isize n_in(dim / 2);
+    isize n_eq(dim / 2);
+    isize n_in(dim / 2);
     std::cout << "dim: " << dim << " n_eq: " << n_eq << " n_in: " << n_in
               << " box: " << dim << std::endl;
     T strong_convexity_factor(1.e-2);
 
-    proxqp::dense::Model<T> qp_random = proxqp::utils::dense_strongly_convex_qp(
+    common::dense::Model<T> qp_random = common::utils::dense_strongly_convex_qp(
       dim, n_eq, n_in, sparsity_factor, strong_convexity_factor);
     Eigen::Matrix<T, Eigen::Dynamic, 1> x_sol =
-      utils::rand::vector_rand<T>(dim);
+      common::utils::rand::vector_rand<T>(dim);
     Eigen::Matrix<T, Eigen::Dynamic, 1> delta(n_in);
-    for (proxqp::isize i = 0; i < n_in; ++i) {
-      delta(i) = utils::rand::uniform_rand();
+    for (isize i = 0; i < n_in; ++i) {
+      delta(i) = common::utils::rand::uniform_rand();
     }
     qp_random.u = qp_random.C * x_sol + delta;
     qp_random.b = qp_random.A * x_sol;
     Eigen::Matrix<T, Eigen::Dynamic, 1> u_box(dim);
     u_box.setZero();
     Eigen::Matrix<T, Eigen::Dynamic, 1> l_box(dim);
     l_box.setZero();
-    for (proxqp::isize i = 0; i < dim; ++i) {
-      T shift = utils::rand::uniform_rand();
+    for (isize i = 0; i < dim; ++i) {
+      T shift = common::utils::rand::uniform_rand();
       u_box(i) = x_sol(i) + shift;
       l_box(i) = x_sol(i) - shift;
     }

benchmark/timings-dense-backend.cpp

@@ -3,13 +3,13 @@
 //
 #include <iostream>
 #include <proxsuite/proxqp/dense/dense.hpp>
-#include <proxsuite/proxqp/utils/random_qp_problems.hpp>
+#include <proxsuite/common/utils/random_qp_problems.hpp>
 
 using T = double;
 using I = long long;
 
 using namespace proxsuite;
-using namespace proxsuite::proxqp;
+using namespace proxsuite::common;
 
 int
 main(int /*argc*/, const char** /*argv*/)
@@ -20,32 +20,32 @@ main(int /*argc*/, const char** /*argv*/)
   T sparsity_factor = 0.75;
   T eps_abs = T(1e-9);
   T elapsed_time = 0.0;
-  proxqp::utils::rand::set_seed(1);
+  common::utils::rand::set_seed(1);
   std::cout << "Dense QP" << std::endl;
-  for (proxqp::isize dim = 100; dim <= 1000; dim = dim + 100) {
+  for (isize dim = 100; dim <= 1000; dim = dim + 100) {
 
-    proxqp::isize n_eq(dim * 2);
-    proxqp::isize n_in(dim * 2);
+    isize n_eq(dim * 2);
+    isize n_in(dim * 2);
     std::cout << "dim: " << dim << " n_eq: " << n_eq << " n_in: " << n_in
               << " box: " << dim << std::endl;
     T strong_convexity_factor(1.e-2);
 
-    proxqp::dense::Model<T> qp_random = proxqp::utils::dense_strongly_convex_qp(
+    common::dense::Model<T> qp_random = common::utils::dense_strongly_convex_qp(
       dim, n_eq, n_in, sparsity_factor, strong_convexity_factor);
     Eigen::Matrix<T, Eigen::Dynamic, 1> x_sol =
-      utils::rand::vector_rand<T>(dim);
+      common::utils::rand::vector_rand<T>(dim);
     Eigen::Matrix<T, Eigen::Dynamic, 1> delta(n_in);
-    for (proxqp::isize i = 0; i < n_in; ++i) {
-      delta(i) = utils::rand::uniform_rand();
+    for (isize i = 0; i < n_in; ++i) {
+      delta(i) = common::utils::rand::uniform_rand();
     }
     qp_random.u = qp_random.C * x_sol + delta;
     qp_random.b = qp_random.A * x_sol;
     Eigen::Matrix<T, Eigen::Dynamic, 1> u_box(dim);
     u_box.setZero();
     Eigen::Matrix<T, Eigen::Dynamic, 1> l_box(dim);
     l_box.setZero();
-    for (proxqp::isize i = 0; i < dim; ++i) {
-      T shift = utils::rand::uniform_rand();
+    for (isize i = 0; i < dim; ++i) {
+      T shift = common::utils::rand::uniform_rand();
       u_box(i) = x_sol(i) + shift;
       l_box(i) = x_sol(i) - shift;
     }
@@ -64,7 +64,9 @@ main(int /*argc*/, const char** /*argv*/)
 
     elapsed_time = 0.0;
     timer.stop();
-    proxqp::dense::QP<T> qp{ dim, n_eq, n_in, true, DenseBackend::PrimalLDLT };
+    proxqp::dense::QP<T> qp{
+      dim, n_eq, n_in, true, common::DenseBackend::PrimalLDLT
+    };
     qp.settings.eps_abs = eps_abs;
     qp.settings.eps_rel = 0;
     // qp.settings.verbose = true;
@@ -97,7 +99,7 @@ main(int /*argc*/, const char** /*argv*/)
 
     elapsed_time = 0.0;
     proxqp::dense::QP<T> qp_compare{
-      dim, n_eq, n_in, true, DenseBackend::PrimalDualLDLT
+      dim, n_eq, n_in, true, common::DenseBackend::PrimalDualLDLT
     };
     qp_compare.settings.eps_abs = eps_abs;
     qp_compare.settings.eps_rel = 0;
@@ -131,7 +133,7 @@ main(int /*argc*/, const char** /*argv*/)
               << elapsed_time * 1e-3 / smooth << "ms" << std::endl;
     elapsed_time = 0.0;
     proxqp::dense::QP<T> qp_compare_bis{
-      dim, n_eq, n_in, true, DenseBackend::Automatic
+      dim, n_eq, n_in, true, common::DenseBackend::Automatic
     };
     qp_compare_bis.settings.eps_abs = eps_abs;
     qp_compare_bis.settings.eps_rel = 0;

benchmark/timings-diagonal-hessian.cpp

@@ -3,13 +3,13 @@
 //
 #include <iostream>
 #include <proxsuite/proxqp/dense/dense.hpp>
-#include <proxsuite/proxqp/utils/random_qp_problems.hpp>
+#include <proxsuite/common/utils/random_qp_problems.hpp>
 
 using T = double;
 using I = long long;
 
 using namespace proxsuite;
-using namespace proxsuite::proxqp;
+using namespace proxsuite::common;
 
 int
 main(int /*argc*/, const char** /*argv*/)
@@ -20,32 +20,32 @@ main(int /*argc*/, const char** /*argv*/)
   T sparsity_factor = 0.75;
   T eps_abs = T(1e-9);
   T elapsed_time = 0.0;
-  proxqp::utils::rand::set_seed(1);
+  common::utils::rand::set_seed(1);
   std::cout << "Dense QP" << std::endl;
-  for (proxqp::isize dim = 100; dim <= 500; dim = dim + 100) {
+  for (isize dim = 100; dim <= 500; dim = dim + 100) {
 
-    proxqp::isize n_eq(dim / 2);
-    proxqp::isize n_in(dim / 2);
+    isize n_eq(dim / 2);
+    isize n_in(dim / 2);
     std::cout << "dim: " << dim << " n_eq: " << n_eq << " n_in: " << n_in
               << " box: " << dim << std::endl;
     T strong_convexity_factor(1.e-2);
 
-    proxqp::dense::Model<T> qp_random = proxqp::utils::dense_strongly_convex_qp(
+    common::dense::Model<T> qp_random = common::utils::dense_strongly_convex_qp(
       dim, n_eq, n_in, sparsity_factor, strong_convexity_factor);
     Eigen::Matrix<T, Eigen::Dynamic, 1> x_sol =
-      utils::rand::vector_rand<T>(dim);
+      common::utils::rand::vector_rand<T>(dim);
     Eigen::Matrix<T, Eigen::Dynamic, 1> delta(n_in);
-    for (proxqp::isize i = 0; i < n_in; ++i) {
-      delta(i) = utils::rand::uniform_rand();
+    for (isize i = 0; i < n_in; ++i) {
+      delta(i) = common::utils::rand::uniform_rand();
     }
     qp_random.u = qp_random.C * x_sol + delta;
     qp_random.b = qp_random.A * x_sol;
     Eigen::Matrix<T, Eigen::Dynamic, 1> u_box(dim);
     u_box.setZero();
     Eigen::Matrix<T, Eigen::Dynamic, 1> l_box(dim);
     l_box.setZero();
-    for (proxqp::isize i = 0; i < dim; ++i) {
-      T shift = utils::rand::uniform_rand();
+    for (isize i = 0; i < dim; ++i) {
+      T shift = common::utils::rand::uniform_rand();
       u_box(i) = x_sol(i) + shift;
       l_box(i) = x_sol(i) - shift;
     }
@@ -69,12 +69,9 @@ main(int /*argc*/, const char** /*argv*/)
 
     elapsed_time = 0.0;
     timer.stop();
-    proxqp::dense::QP<T> qp{ dim,
-                             n_eq,
-                             n_in,
-                             true,
-                             proxsuite::proxqp::DenseBackend::PrimalDualLDLT,
-                             proxsuite::proxqp::HessianType::Diagonal };
+    proxqp::dense::QP<T> qp{
+      dim, n_eq, n_in, true, DenseBackend::PrimalDualLDLT, HessianType::Diagonal
+    };
     qp.settings.eps_abs = eps_abs;
     qp.settings.eps_rel = 0;
     // qp.settings.verbose = true;
@@ -107,12 +104,7 @@ main(int /*argc*/, const char** /*argv*/)
 
     elapsed_time = 0.0;
     proxqp::dense::QP<T> qp_compare{
-      dim,
-      n_eq,
-      n_in,
-      true,
-      proxsuite::proxqp::DenseBackend::PrimalDualLDLT,
-      proxsuite::proxqp::HessianType::Dense
+      dim, n_eq, n_in, true, DenseBackend::PrimalDualLDLT, HessianType::Dense
     };
     qp_compare.settings.eps_abs = eps_abs;
     qp_compare.settings.eps_rel = 0;

benchmark/timings-lp.cpp

@@ -3,13 +3,13 @@
 //
 #include <iostream>
 #include <proxsuite/proxqp/dense/dense.hpp>
-#include <proxsuite/proxqp/utils/random_qp_problems.hpp>
+#include <proxsuite/common/utils/random_qp_problems.hpp>
 
 using T = double;
 using I = long long;
 
 using namespace proxsuite;
-using namespace proxsuite::proxqp;
+using namespace proxsuite::common;
 
 int
 main(int /*argc*/, const char** /*argv*/)
@@ -20,33 +20,32 @@ main(int /*argc*/, const char** /*argv*/)
   T sparsity_factor = 0.75;
   T eps_abs = T(1e-9);
   T elapsed_time = 0.0;
-  proxqp::utils::rand::set_seed(1);
+  common::utils::rand::set_seed(1);
   std::cout << "Dense QP" << std::endl;
-  for (proxqp::isize dim = 10; dim <= 1000;
-       dim = (dim == 10) ? 100 : dim + 100) {
+  for (isize dim = 10; dim <= 1000; dim = (dim == 10) ? 100 : dim + 100) {
 
     if (dim == 10 || dim == 100) {
       smooth = 1000;
     } else {
       smooth = 100;
     }
 
-    proxqp::isize n_eq(dim / 2);
-    proxqp::isize n_in(dim / 2);
+    isize n_eq(dim / 2);
+    isize n_in(dim / 2);
     T strong_convexity_factor(1.e-2);
     std::cout << "dim: " << dim << " n_eq: " << n_eq << " n_in: " << n_in
               << std::endl;
 
-    proxqp::dense::Model<T> qp_random = proxqp::utils::dense_strongly_convex_qp(
+    common::dense::Model<T> qp_random = common::utils::dense_strongly_convex_qp(
       dim, n_eq, n_in, sparsity_factor, strong_convexity_factor);
     qp_random.H.setZero();
-    auto y_sol = proxqp::utils::rand::vector_rand<T>(n_eq);
+    auto y_sol = common::utils::rand::vector_rand<T>(n_eq);
     qp_random.g = -qp_random.A.transpose() * y_sol;
 
     elapsed_time = 0.0;
     timer.stop();
     proxqp::dense::QP<T> qp{
-      dim, n_eq, n_in, false, proxqp::HessianType::Zero
+      dim, n_eq, n_in, false, common::HessianType::Zero
     };
     qp.settings.eps_abs = eps_abs;
     qp.settings.eps_rel = 0;
@@ -77,7 +76,7 @@ main(int /*argc*/, const char** /*argv*/)
 
     elapsed_time = 0.0;
     proxqp::dense::QP<T> qp_compare{
-      dim, n_eq, n_in, false, proxqp::HessianType::Dense
+      dim, n_eq, n_in, false, common::HessianType::Dense
     };
     qp_compare.settings.eps_abs = eps_abs;
     qp_compare.settings.eps_rel = 0;