GLPK (python swiglpk) "Problem has no primal feasible solution" but ok with CVXPY

Question

I'm trying to solve a simple optimization problem:

max   x+y
s.t.  -x <= -1
       x,y in {0,1}^2

using following code

import swiglpk
import numpy as np

def solve_boolean_lp_swig(obj: np.ndarray, aub: np.ndarray, bub: np.ndarray, minimize: bool) -> tuple:

    """
        Solves following optimization problem
        min/max     obj.dot(x)
        s.t         aub.dot(x) <= bub
                    x \in {0, 1}

        obj : m vector
        aub : nxm matrix
        bub : n vector
    """

    # init problem
    ia = swiglpk.intArray(1+aub.size); ja = swiglpk.intArray(1+aub.size)
    ar = swiglpk.doubleArray(1+aub.size)
    lp = swiglpk.glp_create_prob()

    # set obj to minimize if minimize==True else maximize
    swiglpk.glp_set_obj_dir(lp, swiglpk.GLP_MIN if minimize else swiglpk.GLP_MAX)

    # number of rows and columns as n, m 
    swiglpk.glp_add_rows(lp, int(aub.shape[0]))
    swiglpk.glp_add_cols(lp, int(aub.shape[1]))
    
    # setting row constraints (-inf < x <= bub[i])
    for i, v in enumerate(bub):
        swiglpk.glp_set_row_bnds(lp, i+1, swiglpk.GLP_UP, 0.0, float(v))
    
    # setting column constraints (x in {0, 1})
    for i in range(aub.shape[1]):
        # not sure if this is needed but perhaps for presolving
        swiglpk.glp_set_col_bnds(lp, i+1, swiglpk.GLP_FR, 0.0, 0.0) 
        # setting x in {0,1}
        swiglpk.glp_set_col_kind(lp, i+1, swiglpk.GLP_BV)

    # setting aub 
    for r, (i,j) in enumerate(np.argwhere(aub != 0)):
        ia[r+1] = int(i)+1; ja[r+1] = int(j)+1; ar[r+1] = float(aub[i,j])

    # solver settings
    iocp = swiglpk.glp_iocp()
    swiglpk.glp_init_iocp(iocp)
    iocp.msg_lev = swiglpk.GLP_MSG_ALL
    iocp.presolve = swiglpk.GLP_ON
    iocp.binarize = swiglpk.GLP_ON

    # setting objective
    for i,v in enumerate(obj):
        swiglpk.glp_set_obj_coef(lp, i+1, float(v))
    
    swiglpk.glp_load_matrix(lp, r, ia, ja, ar)
    info = swiglpk.glp_intopt(lp, iocp)
    
    # use later
    #status = swiglpk.glp_mip_status(lp)
    
    x = np.array([swiglpk.glp_mip_col_val(lp, int(i+1)) for i in range(obj.shape[0])])

    # for now, keep it simple. info == 0 means optimal 
    # solution (there are others telling feasible solution)
    return (info == 0), x

and the following instance (as given on top)

solve_boolean_lp_swig(
    obj = np.array([ 1, 1]),
    aub = np.array([[-1, 0]]),
    bub = np.array([-1]),
    minimize = False
)

In my mind x=[1,0] should be a valid solution since dot([-1, 0], x) <= -1 (and [1,0] are boolean) holds but solver says PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION. However, if i run the same problem instance using the lib CVXOPT instead, with cvxopt.glpk.ilp, the solver finds an optimal solution. I've seen the c-code underneath cvxopt and has done the same so I suspect something small that I cannot see..

Answer 1

Add to the model:

swiglpk.glp_write_lp(lp,None,"xxx.lp")

Then you'll see immediately what the problem is:

\* Problem: Unknown *\

Maximize
 obj: + z_1 + z_2

Subject To
 r_1: 0 z_1 <= -1

Bounds
 0 <= z_1 <= 1
 0 <= z_2 <= 1

Generals
 z_1
 z_2

End

I noticed that r=0, so the ne argument for the load call is already wrong. If you set r=1 things look better.