Possible Bug in Setting Contraints

[reichardtj]

First off, I'd like to say thanks to Gianluca and all contributors for supplying this really amazing piece of software. Great effort!

In playing with the python interface, I noticed that hpipm.d_ocp_qp_set("idxbx", np.arange(number_of_states), 0) gives me an error while hpipm.d_ocp_qp_set("Jbx", np.eye(number_of_states), 0) does not and my trials work.

Digging deeper, I found the following lines that all feature an odd:

if (nbx[ii>0])

hpipm/ocp_qp/x_ocp_qp.c

Line 556 in f257970

if(nbx[ii>0])

hpipm/ocp_qp/x_ocp_qp.c

Line 563 in f257970

if(nb[ii>0])

hpipm/ocp_qp/x_ocp_qp.c

Line 653 in f257970

if(nbx[ii>0])

hpipm/ocp_qp/x_ocp_qp.c

Line 660 in f257970

if(nb[ii>0])

I admit that I am not fully clear on the logic of the code and the memory layout, but I suspect it should read

if (nbx[ii] > 0)

Thanks for considering this issue.
Best
Jörg

[giaf]

First of all, thanks for appreciating HPIPM :)
Yes you are right with the bug report, thanks for spotting and finding the exact issue, I'll fix now!

[giaf]

The issue you report should be fixed here db2fe4d

However, this should not affect the function d_ocp_qp_set() that you mention was failing for the field idxbx...
I guess there is still something else...

[reichardtj]

Excellent! Thanks for the quick fix. It does not solve the problem with set('idxbx', ...) though as you already expected. I've reinstalled after your fix from both the latest master of Blasfeo and HPIPM, however, the problem persists.

Here is a minimal example showing the problem. From looking at the python interface, I doubt it's a python issue:

import os
# this is needed only on windows where the library searches the path instead of LD_LIBRARY_PATH
#os.environ['PATH'] = os.environ['PATH'] + 'D:\\hpipm\\lib;D:\\github_repos\\blasfeo\\lib'

# On Mac and presumably on Linux this is needed instead to find the library
from ctypes import CDLL
_ = CDLL('/opt/hpipm/lib/libhpipm.so')

import numpy as np

import hpipm_python as hpipm

T = 5

u_lower = np.array([-2, -1])
u_upper = np.array([ 2,  1])

A = np.array([
    [1, 0, 0.1, 0],
    [0, 1, 0, 0.1],
    [0, 0, 1, 0],
    [0, 0, 0, 1]
])
B = np.array([
    [0, 0],
    [0, 0],
    [0.1, 0],
    [0, 0.1]
])

nx = A.shape[0]
nu = B.shape[1]

dim = hpipm.hpipm_ocp_qp_dim(T)

dim.set('nx', nx, 0, T)   # number of states
dim.set('nu', nu, 0, T-1) # number of inputs

dim.set('nbx', 0, 0, T)
dim.set('nbx', nx, 0)     # initial condition

dim.set('nbu', nu, 0, T-1)  # for upper and lower bounds on u

qp = hpipm.hpipm_ocp_qp(dim)

qp_sol = hpipm.hpipm_ocp_qp_sol(dim)

mode = 'speed'
arg = hpipm.hpipm_ocp_qp_solver_arg(dim, mode)

# set up solver
solver = hpipm.hpipm_ocp_qp_solver(dim, arg)

#qp.set('idxbx', np.arange(nx).astype(int), 0) ## ====> this does not work
qp.set('Jbx', np.eye(nx), 0) ## ======> This works

qp.set('lbx', np.zeros(nx), 0)
qp.set('ubx', np.zeros(nx), 0)

qp.set('idxbu', np.arange(nu).astype(int), 0, T-1) #=====> this does work
#qp.set('Jbu', np.eye(nu), 0, T-1)  ## This works, too!
qp.set('lbu', u_lower, 0, T-1)
qp.set('ubu', u_upper, 0, T-1)

x = np.zeros((T + 1, nx))
u = np.clip(np.random.randn(T, nu) * 0.1, u_lower, u_upper)

x[0] = np.array([1, 1, 0, 0])
for i in range(T):
    x[i + 1] = A @ x[i] + B @ u[i]
    qp.set('A', A, i)
    qp.set('B', B, i)
    qp.set('b', np.zeros(nx), i)

    qp.set('Q', np.zeros((nx, nx)), i)
    qp.set('S', np.zeros((nu, nx)), i)
    qp.set('R', np.eye(nu) * 0.01, i)
    qp.set('r', np.zeros(nu), i)
    qp.set('q', np.zeros(nx), i)

    qp.set('lbu', u_lower - u[i], i)
    qp.set('ubu', u_upper - u[i], i)


qp.set('Q', np.eye(nx), T)
qp.set('S', np.zeros((nu, nx)), T)
qp.set('R', np.zeros((nu, nu)), T)
qp.set('r', np.zeros(nu), T)
qp.set('q', np.zeros(nx), T)


# solve qp    
solver.solve(qp, qp_sol)

On my Windows install, the solver crashes with an error message quoting some address it's not supposed to read from. On my Mac, the whole Jupyter Kernel dies without error message. The curious thing is that set("idxbu",...) does not not cause a problem, only set("idxbx",...).

As a first check, it would be good to see if other people can reproduce the problem.

[giaf]

Hi, thanks for the example, it was very useful to find the bug!

The issue was that the integer arrays were converted into double arrays in the python interface, this should fix it 672150d
BTW it was wrong also for idxbu, since internally the C code was getting [0 0] instead of [0 1] in this case.

[reichardtj]

Great. I've tried that and that all seems to work. Thanks for fixing this so quickly!