hpipm_dense_qp_solver fails with only equality constraints

[steradiant]

Hello, I intend to use HPIPM for solving a quadratic static constrained optimization problem. I use the MATLAB API, more accurate, the hpipm_dense_qp_solver function. My problem looks as follows

min_x 1/2x'Hx s.t. Ax=b

According to the publication https://arxiv.org/abs/2003.02547 I set the following arguments

dim = hpipm_dense_qp_dim();
dim.set('nv', nv); % nv is the number of optimization variables, in my case 94
dim.set('ne', ne); % ne is the number of equality constraints, in my case 89
qp = hpipm_dense_qp(dim);
qp.set('H', H);
qp.set('A', Aeq);
qp.set('b', beq);
sol = hpipm_dense_qp_sol(dim);
mode = 'speed';
arg = hpipm_dense_qp_solver_arg(dim, mode);
    
solver = hpipm_dense_qp_solver(dim, arg);

solver.solve(qp, sol);
x = sol.get('v');

The solver exits with status 0 but the values of the solution are all NaN. If I try to solve the same problem with quadprog from MATLAB, I get a valid solution, hence, the problem is solvable. Where could the problem be?

[giaf]

The issue could be in the problem being badly/ill conditioned, or a bug in the solver. It's difficult to say without having the problem data. You can try to check the rank/conditioning of the H and Aeq matrices. If all looks fine there, and if you can share the data, I could find some time at the end of next week to check it.

[steradiant]

Because it's current research, it's pretty hard to share more detailed data. However, I tried the following formulation:

dim = hpipm_dense_qp_dim();
dim.set('nv', nv);
dim.set('ng', ne);
qp = hpipm_dense_qp(dim);
qp.set('H', H);
qp.set('C', Aeq);
qp.set('lg', beq); 
qp.set('ug', beq);

sol = hpipm_dense_qp_sol(dim);
mode = 'speed';
arg = hpipm_dense_qp_solver_arg(dim, mode);
    
solver = hpipm_dense_qp_solver(dim, arg);
    
solver.solve(qp, sol);
    
x = sol.get('v');

This problem is solved successfully. Therefore I'd guess, that there is an error in the implementation of the equality constraints. Are all variables named like in the publication? I sadly can't find a reference on how the names of the publication relate to the names in MATLAB (or more precisely, in the get and set functions)? Also, an explanation of the terminal output would be great, e.g. the meaning of idxb, d_mask, ...

[alexphieler]

Hi, I think I am facing as similar issue. With the following formulation:

int hpipm_status;

int dim_size = d_dense_qp_dim_memsize();
void *dim_mem = malloc(dim_size);

struct d_dense_qp_dim dim;

d_dense_qp_dim_create(&dim, dim_mem);

d_dense_qp_dim_set_all(6, 2, 0, 0, 0, 0, &dim);


// /************************************************
// * qp
// ************************************************/

int qp_size = d_dense_qp_memsize(&dim);
void *qp_mem = malloc(qp_size);

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

double w = 0.01;
double wA = 10.0;
double wR = 10.0;
double targetA = 5.0;
double targetR = 0.0;
double t2f = 13.95 * 2.0 / 0.47;
double al = 1.0;
double ar = 1.0;
double iZ = 110.0;
double m = 200;
double lF = 0.75;
double tW = 1.15;


double H[] = {
                    w, 0.0, 0.0, 0.0, 0.0, 0.0,
                    0.0, w, 0.0, 0.0, 0.0, 0.0,
                    0.0, 0.0, w, 0.0, 0.0, 0.0,
                    0.0, 0.0, 0.0, w, 0.0, 0.0,
                    0.0, 0.0, 0.0, 0.0, wA, 0.0,
                    0.0, 0.0, 0.0, 0.0, 0.0, wR
                    };

double g[] = {0.0, 0.0, 0.0, 0.0, -wA*targetA, -wR*targetR};

double A[] = {
                    t2f*1, t2f*1, t2f*al, t2f*ar, -1*m, 0.0,
                    -0.5*tW*t2f*1, 0.5*tW*t2f*1, -0.5*tW*t2f*al, 0.5*tW*t2f*al, 0.0, -1*iZ
			};

double b[] = {0.0, 0.0};


int idxb[] = {};
double d_lb[] = {};
double d_ub[] = {};

double C[] = {};
double d_lg[] = {};
double d_ug[] = {};

double Zl[] = {};
double Zu[] = {};

double zl[] = {};
double zu[] = {};

int idxs[] = {};

double d_ls[] = {};
double d_us[] = {};

d_dense_qp_set_all(H, g, A, b, idxb, d_lb, d_ub, C, d_lg, d_ug, Zl, Zu, zl, zu, idxs, d_ls, d_us, &qp);


// /************************************************
// * qp sol
// ************************************************/

int qp_sol_size = d_dense_qp_sol_memsize(&dim);
void *qp_sol_mem = malloc(qp_sol_size);

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

// /************************************************
// * ipm arg
// ************************************************/

int ipm_arg_size = d_dense_qp_ipm_arg_memsize(&dim);
void *ipm_arg_mem = malloc(ipm_arg_size);

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

d_dense_qp_ipm_arg_set_default(mode, &arg);


// /************************************************
// * ipm workspace
// ************************************************/

int ipm_size = d_dense_qp_ipm_ws_memsize(&dim, &arg);
void *ipm_mem = malloc(ipm_size);

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

// /************************************************
// * ipm solver
// ************************************************/

d_dense_qp_ipm_solve(&qp, &qp_sol, &arg, &workspace);

The solver returns with status 0, while the vector v of the solution is: -1.11906 -1.11906 -1.33013 -1.09747 4.99890 0.00281 which is definitely not correct.

As comparison, when I run the same problem with OSQP, I get [4.20550227 4.20550227 4.20550227 4.20550227 4.99291546 0. ]

A reformulation of the equality constraint as an inequality constraint with equal upper and lower limits did also not solve the issue for me. The problem occurs on both the "master" and "stable" branch.

[giaf]

Hi, maybe there is an issue with the problem setup. If I print out the matrices and vectors in the internal QP structure (using d_dense_qp_print) and use e.g. octave to solve the same QP, I get the same solution as HPIPM:

octave:7> H
H =

    0.0100         0         0         0         0         0
         0    0.0100         0         0         0         0
         0         0    0.0100         0         0         0
         0         0         0    0.0100         0         0
         0         0         0         0   10.0000         0
         0         0         0         0         0   10.0000

octave:8> q
q =

    0
    0
    0
    0
  -50
    0

octave:9> A
A =

    59.3617    59.3617  -200.0000   -34.1330   -34.1330          0
    59.3617    59.3617          0    34.1330    34.1330  -110.0000

octave:10> b
b =

   0
   0

octave:11> qp([], H, q, A, b)
ans =

  -1.1191e+00
  -1.1191e+00
  -1.3301e+00
  -1.0975e+00
   4.9989e+00
   2.8052e-03

Could you double-check that this is the QP that you intend to solve?

[giaf]

Maybe an hint: double-check that A is what you intend it to be, and keep in mind that the matrices passed to HPIPM are expected to be in column-major by default.

[alexphieler]

Hi, my problem was indeed caused by the fact that haven't the used column-major format, now it works as expected. Thanks very much for your help!