giaf / hpipm

High-performance interior-point-method QP and QCQP solvers
Other
554 stars 130 forks source link

Equality Constraints not working with dense qp #144

Open dlcfort opened 1 year ago

dlcfort commented 1 year ago

The optimization that I set up does not work when I set the number of general inequality constraints "ng" to zero. However when I set "ng" to a value of 1 it works fine. Is there something wrong with my formulation or is there a bug in the code? You can see from the printed vector (Av) at the end that it does not match up with the vector (b)

from hpipm_python import *
from hpipm_python.common import *
import numpy as np
import time

nv = 8 #number of variables
ne = 8 # number of equality constraints
ng = 0 
dim = hpipm_dense_qp_dim()
dim.set('nv', nv)
dim.set('ne', ne)
dim.set('ng', ng)
H = np.array([[ 2, -2,  0,  0,  0,  0,  0,  0],
                [-2,  4, -2,  0,  0,  0,  0,  0],
                [ 0, -2,  4, -2,  0,  0,  0,  0],
                [ 0,  0, -2,  2,  0,  0,  0,  0],
                [ 0,  0,  0,  0,  2, -2,  0,  0],
                [ 0,  0,  0,  0, -2,  4, -2,  0],
                [ 0,  0,  0,  0,  0, -2,  4, -2],
                [ 0,  0,  0,  0,  0,  0, -2,  2]])
g = np.array([[0.],[0.],[0.],[0.],[0.],[0.],[0.],[0.]])
A = np.array([[ 1,  0,  0,  0,  0,  0,  0,  0],
    [ 0,  0,  0,  0,  1,  0,  0,  0],
    [ 0,  0,  0,  1,  0,  0,  0,  0],
    [ 0,  0,  0,  0,  0,  0,  0,  1],
    [-3,  3,  0,  0,  0,  0,  0,  0],
    [ 0,  0,  0,  0, -3,  3,  0,  0],
    [ 0,  0, -3,  3,  0,  0,  0,  0],
    [ 0,  0,  0,  0,  0,  0, -3,  3]])
b = np.array([[ 1. ],[ 2. ],[10. ],[12. ],[ 1.5],[ 3. ],[ 0. ],[ 1.5]])
qp = hpipm_dense_qp(dim)
qp.set('H', H)
qp.set('g', g)
qp.set('A', A)
qp.set('b', b)
qp_sol = hpipm_dense_qp_sol(dim)
mode = 'robust'
arg = hpipm_dense_qp_solver_arg(dim, mode)
arg.set('mu0', 1e4)
arg.set('iter_max', 30)
arg.set('tol_stat', 1e-4)
arg.set('tol_eq', 1e-5)
arg.set('tol_ineq', 1e-5)
arg.set('tol_comp', 1e-5)
arg.set('reg_prim', 1e-12)
solver = hpipm_dense_qp_solver(dim, arg)
optimization_start_time = time.time()
solver.solve(qp, qp_sol)
optimization_end_time = time.time()
v = qp_sol.get('v')

print("A @ v[:,None]: " , np.dot(A , v))

A @ v[:,None]: [[-6.32751741] [-9.52767775] [40.07520835] [40.07520835] [ 1.5 ] [ 3. ] [68.85408864] [72.90432914]]

giaf commented 1 year ago

Hi, it seems that the issue is caused by the fact that H is rank deficient (two eigenvalues are zero) so the problem is not strictly convex, if you put H to something strictly positive definite then it works fine.

In case of putting ng>0, this creates an inequality constraint with value 0 that is always satisfied, but it triggers the IPM machinery that using some regularization overcomes the non-strictly-convex issue. If there are only equality constraints without inequality, the IPM machinery is not used but only a linear system solve, that apparently fails if H is not strictly positive definite (and thus not invertible too).