|
| 1 | +// |
| 2 | +// Copyright (c) 2023 INRIA |
| 3 | +// |
| 4 | +#include <iostream> |
| 5 | +#include <proxsuite/proxqp/dense/dense.hpp> |
| 6 | +#include <proxsuite/proxqp/utils/random_qp_problems.hpp> |
| 7 | + |
| 8 | +using T = double; |
| 9 | +using I = long long; |
| 10 | + |
| 11 | +using namespace proxsuite; |
| 12 | +using namespace proxsuite::proxqp; |
| 13 | + |
| 14 | +int |
| 15 | +main(int /*argc*/, const char** /*argv*/) |
| 16 | +{ |
| 17 | + Timer<T> timer; |
| 18 | + int smooth = 100; |
| 19 | + |
| 20 | + T sparsity_factor = 0.75; |
| 21 | + T eps_abs = T(1e-9); |
| 22 | + T elapsed_time = 0.0; |
| 23 | + proxqp::utils::rand::set_seed(1); |
| 24 | + std::cout << "Dense QP" << std::endl; |
| 25 | + for (proxqp::isize dim = 100; dim <= 1000; dim = dim + 100) { |
| 26 | + |
| 27 | + proxqp::isize n_eq(dim / 2); |
| 28 | + proxqp::isize n_in(dim / 2); |
| 29 | + std::cout << "dim: " << dim << " n_eq: " << n_eq << " n_in: " << n_in |
| 30 | + << " box: " << dim << std::endl; |
| 31 | + T strong_convexity_factor(1.e-2); |
| 32 | + |
| 33 | + proxqp::dense::Model<T> qp_random = proxqp::utils::dense_strongly_convex_qp( |
| 34 | + dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); |
| 35 | + Eigen::Matrix<T, Eigen::Dynamic, 1> x_sol = |
| 36 | + utils::rand::vector_rand<T>(dim); |
| 37 | + Eigen::Matrix<T, Eigen::Dynamic, 1> delta(n_in); |
| 38 | + for (proxqp::isize i = 0; i < n_in; ++i) { |
| 39 | + delta(i) = utils::rand::uniform_rand(); |
| 40 | + } |
| 41 | + qp_random.u = qp_random.C * x_sol + delta; |
| 42 | + qp_random.b = qp_random.A * x_sol; |
| 43 | + Eigen::Matrix<T, Eigen::Dynamic, 1> u_box(dim); |
| 44 | + u_box.setZero(); |
| 45 | + Eigen::Matrix<T, Eigen::Dynamic, 1> l_box(dim); |
| 46 | + l_box.setZero(); |
| 47 | + for (proxqp::isize i = 0; i < dim; ++i) { |
| 48 | + T shift = utils::rand::uniform_rand(); |
| 49 | + u_box(i) = x_sol(i) + shift; |
| 50 | + l_box(i) = x_sol(i) - shift; |
| 51 | + } |
| 52 | + using Mat = |
| 53 | + Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>; |
| 54 | + Mat C_enlarged(dim + n_in, dim); |
| 55 | + C_enlarged.setZero(); |
| 56 | + C_enlarged.topLeftCorner(n_in, dim) = qp_random.C; |
| 57 | + C_enlarged.bottomLeftCorner(dim, dim).diagonal().array() += 1.; |
| 58 | + Eigen::Matrix<T, Eigen::Dynamic, 1> u_enlarged(n_in + dim); |
| 59 | + Eigen::Matrix<T, Eigen::Dynamic, 1> l_enlarged(n_in + dim); |
| 60 | + u_enlarged.head(n_in) = qp_random.u; |
| 61 | + u_enlarged.tail(dim) = u_box; |
| 62 | + l_enlarged.head(n_in) = qp_random.l; |
| 63 | + l_enlarged.tail(dim) = l_box; |
| 64 | + |
| 65 | + elapsed_time = 0.0; |
| 66 | + timer.stop(); |
| 67 | + proxqp::dense::QP<T> qp{ dim, n_eq, n_in, true }; |
| 68 | + qp.settings.eps_abs = eps_abs; |
| 69 | + qp.settings.eps_rel = 0; |
| 70 | + qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; |
| 71 | + for (int j = 0; j < smooth; j++) { |
| 72 | + timer.start(); |
| 73 | + qp.init(qp_random.H, |
| 74 | + qp_random.g, |
| 75 | + qp_random.A, |
| 76 | + qp_random.b, |
| 77 | + qp_random.C, |
| 78 | + qp_random.l, |
| 79 | + qp_random.u, |
| 80 | + l_box, |
| 81 | + u_box); |
| 82 | + qp.solve(); |
| 83 | + timer.stop(); |
| 84 | + elapsed_time += timer.elapsed().user; |
| 85 | + if (qp.results.info.pri_res > eps_abs || |
| 86 | + qp.results.info.dua_res > eps_abs) { |
| 87 | + std::cout << "dual residual " << qp.results.info.dua_res |
| 88 | + << "; primal residual " << qp.results.info.pri_res |
| 89 | + << std::endl; |
| 90 | + std::cout << "total number of iteration: " << qp.results.info.iter |
| 91 | + << std::endl; |
| 92 | + } |
| 93 | + } |
| 94 | + std::cout << "timings QP with box constraints feature : \t" |
| 95 | + << elapsed_time * 1e-3 / smooth << " ms" << std::endl; |
| 96 | + |
| 97 | + elapsed_time = 0.0; |
| 98 | + proxqp::dense::QP<T> qp_compare{ dim, n_eq, dim + n_in, false }; |
| 99 | + qp_compare.settings.eps_abs = eps_abs; |
| 100 | + qp_compare.settings.eps_rel = 0; |
| 101 | + qp_compare.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; |
| 102 | + for (int j = 0; j < smooth; j++) { |
| 103 | + timer.start(); |
| 104 | + qp_compare.init(qp_random.H, |
| 105 | + qp_random.g, |
| 106 | + qp_random.A, |
| 107 | + qp_random.b, |
| 108 | + C_enlarged, |
| 109 | + l_enlarged, |
| 110 | + u_enlarged); |
| 111 | + qp_compare.solve(); |
| 112 | + timer.stop(); |
| 113 | + elapsed_time += timer.elapsed().user; |
| 114 | + |
| 115 | + if (qp_compare.results.info.pri_res > eps_abs || |
| 116 | + qp_compare.results.info.dua_res > eps_abs) { |
| 117 | + std::cout << "dual residual " << qp_compare.results.info.dua_res |
| 118 | + << "; primal residual " << qp_compare.results.info.pri_res |
| 119 | + << std::endl; |
| 120 | + std::cout << "total number of iteration: " |
| 121 | + << qp_compare.results.info.iter << std::endl; |
| 122 | + } |
| 123 | + } |
| 124 | + std::cout << "timings QP without box constraints feature : \t" |
| 125 | + << elapsed_time * 1e-3 / smooth << " ms" << std::endl; |
| 126 | + } |
| 127 | +} |
0 commit comments