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..
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.