hpipm_dense_qp_solver fails with only box constraints

[driftingsn0w]

Hello, I intend to use HPIPM to solve an the first example from https://www.mathworks.com/help/optim/ug/quadprog.html . The form is H = [1 -1; -1 2]; f = [-2; -6]; A = [1 1; -1 2; 2 1]; b = [2; 2; 3]; There is no equality constraints and only box constraints. My code is as follows, I can't get correct solution and get return code 2. By the way, I can get correct solution using cvxopt with python. Where could the problem be ? Thanks in advance.

`#include

include

include

include "hpipm_d_dense_qp_dim.h"

include "hpipm_d_dense_qp.h"

include "hpipm_d_dense_qp_ipm.h"

/ prints a matrix in column-major format / void d_print_mat(int m, int n, double A, int lda) { int i, j; for(i=0; i<m; i++) { for(j=0; j<n; j++) { printf("%9.5f ", A[i+ldaj]); } printf("\n"); } printf("\n"); }

int main() {

//qp dim and data

int nv = 2;
int ne = 0;
int nb = 3;
int ng = 0;
int ns = 0;
int nsb = 0;
int nsg = 0;

double H[] = {1.0, -1.0, -1.0, 2.0};
double g[] = {-2.0, -6.0};
double A[] = {};
double b[] = {};
double C[] = {1.0, 1.0, -1.0, 2.0, 2.0, 1.0};
double d_ub[] = {2.0, 2.0, 3.0};
double d_lb[] = {};
int idxb[] = {0,1,2};

// double d_lg[] = {};
// double d_ug[] = {};
// double Zl[] = {1e3, 1e3};
// double Zu[] = {1e3, 1e3};
// double zl[] = {1e2, 1e2};
// double zu[] = {1e2, 1e2};
// int idxs[] = {0, 1};
// double d_ls[] = {0, 0};
// double d_us[] = {0, 0};

// dense qp dim
hpipm_size_t dense_qp_dim_size = d_dense_qp_dim_memsize();
void* qp_dim_mem = malloc(dense_qp_dim_size);

struct d_dense_qp_dim qp_dim;
d_dense_qp_dim_create(&qp_dim,qp_dim_mem);
d_dense_qp_dim_set_all(nv, ne, nb, ng, nsb, nsg, &qp_dim);

// dense qp

hpipm_size_t qp_size = d_dense_qp_memsize(&qp_dim);
void* qp_mem = malloc(qp_size);

struct d_dense_qp qp;
d_dense_qp_create(&qp_dim, &qp, qp_mem);

d_dense_qp_set_H(H, &qp);
d_dense_qp_set_g(g, &qp);
d_dense_qp_set_A(A, &qp);
d_dense_qp_set_b(b, &qp);
d_dense_qp_set_idxb(idxb, &qp);
d_dense_qp_set_lb(d_lb, &qp);
d_dense_qp_set_ub(d_ub, &qp);
d_dense_qp_set_C(C, &qp);
// d_dense_qp_set_lg(d_lg, &qp);
// d_dense_qp_set_ug(d_ug, &qp);
// d_dense_qp_set_Zl(Zl, &qp);
// d_dense_qp_set_Zu(Zu, &qp);
// d_dense_qp_set_zl(zl, &qp);
// d_dense_qp_set_zu(zu, &qp);
// d_dense_qp_set_idxs(idxs, &qp);
// d_dense_qp_set_ls(d_ls, &qp);
// d_dense_qp_set_us(d_us, &qp);

// dense qp sol
hpipm_size_t qp_sol_size = d_dense_qp_sol_memsize(&qp_dim);
void* qp_sol_mem = malloc(qp_sol_size);

struct d_dense_qp_sol qp_sol;
d_dense_qp_sol_create(&qp_dim, &qp_sol, qp_sol_mem);

// ipm arg
hpipm_size_t ipm_arg_size = d_dense_qp_ipm_arg_memsize(&qp_dim);
void* ipm_arg_mem = malloc(ipm_arg_size);

struct d_dense_qp_ipm_arg arg;
d_dense_qp_ipm_arg_create(&qp_dim, &arg, ipm_arg_mem);

enum hpipm_mode mode = ROBUST;

d_dense_qp_ipm_arg_set_default(mode, &arg);

int iter_max = 20; //25;
double mu0 = 1e1;
int comp_res_exit = 1;
double tol_stat = 1e-12;
double tol_eq = 1e-12;
double tol_ineq = 1e-12;
double tol_comp = 1e-12;
int kkt_fact_alg = 0;
int remove_lin_dep_eq = 1;
double tau_min = 1e-3;
int test = 1e-5;

d_dense_qp_ipm_arg_set_iter_max(&iter_max, &arg);
d_dense_qp_ipm_arg_set_mu0(&mu0, &arg);
d_dense_qp_ipm_arg_set_comp_res_exit(&comp_res_exit, &arg);
d_dense_qp_ipm_arg_set_tol_stat(&tol_stat, &arg);
d_dense_qp_ipm_arg_set_tol_eq(&tol_eq, &arg);
d_dense_qp_ipm_arg_set_tol_ineq(&tol_ineq, &arg);
d_dense_qp_ipm_arg_set_tol_comp(&tol_comp, &arg);
d_dense_qp_ipm_arg_set_kkt_fact_alg(&kkt_fact_alg, &arg);
d_dense_qp_ipm_arg_set_remove_lin_dep_eq(&remove_lin_dep_eq, &arg);
// ipm
hpipm_size_t ipm_size = d_dense_qp_ipm_ws_memsize(&qp_dim, &arg);

void* ipm_mem = malloc(ipm_size);

struct d_dense_qp_ipm_ws workspace;
d_dense_qp_ipm_ws_create(&qp_dim, &arg, &workspace, ipm_mem);

// d_dense_qp_remove_lin_dep_eq(&qp, &arg, &workspace);

// d_dense_qp_restore_lin_dep_eq(&qp, &arg, &workspace);

int rep, nrep = 1;
int hpipm_return; // 0 normal, 1 max iter

for (rep = 0; rep < nrep;rep++) {
  d_dense_qp_ipm_solve(&qp, &qp_sol, &arg, &workspace);
  d_dense_qp_ipm_get_status(&workspace, &hpipm_return);
}
printf("Dense qp solver return %d \n", hpipm_return);

printf("\nsolution\n\n");
printf("\nv\n");
d_print_mat(1, nv, qp_sol.v->pa, 1);

// free memory

free(qp_dim_mem);
free(qp_mem);
free(qp_sol_mem);
free(ipm_arg_mem);
free(ipm_mem);
// hint
printf("Hpipm successful.\n");

return 0; 

}

`

Output is: Dense qp solver return 2 solution v -0.77934 0.52175

[giaf]

Hi, thanks for trying out HPIPM.

There seem to be a number of issues in your code:

1. The constraints in the example are general constraints, not box constraints, so it should be nb=0 and ng=3. Consequently, the bounds should be passed in d_lg and d_ug instead of d_lb and d_ub; additionally, no idxb should be passed.
2. In the C interface, the data is expected to be passed in column-major, this means that you are actually passing the transposed of the C matrix you intend to. You intended C matrix in column major reads {1.0, -1.0, 2.0, 1.0, 2.0, 1.0}.
3. If no lower bound is passed, then it is assumed to be 0, that is not what you expect in this problem. You should instead pass something like {-1e6, -1e6, -1e6} as lower bound (i.e. l_ug).
4. Standards IPM arguments are totally fine to successfully solve this problem (so no need to switch to ROBUST mode and overwrite standard arguments). Actually with the arguments you choose, the IPM returns status 1 as the exit tolerances set to 1e-12 are too tight. Standard ones or e.g. 1e-8 would make the job.
5. There is a number of compilations errors to fix too. And it is strongly suggested to not access BLASFEO structures directly, but rather using dedicated functions. E.g to print a vector you can use blasfeo_print_dvec https://github.com/giaf/blasfeo/blob/master/include/blasfeo_d_aux_ext_dep.h#L121

Dense qp solver return 0 

solution

v
  0.66667   1.33333 

Hpipm successful.

[driftingsn0w]

Thanks for your help :) My understanding of how to use HPIPM is not accurate enough. According to your suggestion, I have successfully solved the problem.