Failure to solve QP with equality at the boundary

[adamheins]

The dense QP solver (used here via the Python interface) fails on the following example. Notice that the QP is box-constrained with a single equality constraint that requires the first decision variable to be exactly at the lower bound of the box constraints.

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

# dim
nv = 3  # number of variables
ne = 1  # number of equality constraints
nb = 3  # number of box constraints
ng = 0  # number of general (inequality) constraints

dim = hpipm_dense_qp_dim()
dim.set('nv', nv)
dim.set('nb', nb)
dim.set('ne', ne)
dim.set('ng', ng)

# data
H = np.array([[ 65., -22., -16.], [-22., 14., 7.], [-16., 7., 5.]])
g = np.array([-13.,  15.,   7.])
A = np.array([[1., 0, 0]])
b = np.array([-0.5])
lb = np.array([-0.5, -2, -0.8])

# qp
qp = hpipm_dense_qp(dim)
qp.set('H', H)
qp.set('g', g)
qp.set('A', A)
qp.set('b', b)
qp.set('idxb', np.array([0, 1, 2]))
qp.set('lb', lb)
qp.set('ub_mask', np.zeros(3))

arg = hpipm_dense_qp_solver_arg(dim, 'robust')
qp_sol = hpipm_dense_qp_sol(dim)
solver = hpipm_dense_qp_solver(dim, arg)
solver.solve(qp, qp_sol)

The expected result is that the QP is solved with primal solution [0.4, -0.4, 1.0] (this problem is taken from here). However, the actual result I get is

v      = [-0.5        -1.45714286 -0.8       ]
pi     = [-5.68933053e+15]
lam_lb = [5.68933053e+15 1.84210526e-16 8.00000000e-01]
lam_ub = [0. 0. 0.]

solver statistics:
ipm return = 1
ipm max res stat = 1.842105e-16
ipm max res eq   = 0.000000e+00
ipm max res ineq = 2.220446e-16
ipm max res comp = 5.689331e-01
ipm iter = 100

Interestingly, and perhaps related to this comment, converting the single equality constraint Av = b to the double-sided inequality b <= Av <= b does claim success (though the primal optimal point is the same), with result below:

v      = [-0.5        -1.45714286 -0.8       ]
pi     = []
lam_lb = [3.17618472e+02 3.42291874e-13 8.00000000e-01]
lam_ub = [0. 0. 0.]
lam_lg = [317.61847224]
lam_ug = [635.87980162]

solver statistics:
ipm return = 0
ipm max res stat = 2.501110e-12
ipm max res eq   = 0.000000e+00
ipm max res ineq = 1.110223e-16
ipm max res comp = 3.235833e-12
ipm iter = 7

Let me know if you have any insights on this.

[giaf]

Hi, I think you may be mixing problem formulation and solution from different instances.
The solution you propose can't be correct, simply because using the general constraint you constraint x[0]=-0.5 (starting the indexing from 0), but your proposed solution has x[0]=0.4.

Double checking with the qp solver in octave I get the same primal solution as the one returned by HPIPM.

In any case, clearly HPIPM seems to have some numerical issue in the dual solution computation with this instance, likely due to the fact that the equality constraint and the lower bound on x[0] are the same.
I will look into this.

[adamheins]

Apologies, you are quite right: I mixed up the problem I was looking at. The primal solution should indeed be [-0.5, -1.45714286, -0.8], which is the solution that HPIPM does achieve. It just has some numerical difficulties with the dual, as you say. Thanks for taking a look!