How can I get integer solutions with scipy.optimize.linprog?

Solution 1:

From the docs:

method : str, optional Type of solver. At this time only ‘simplex’ is supported.

Simplex cannot handle integrality constraints so you cannot solve integer programming problems with scipy.optimize.linprog yet. You can try other libraries like PuLP, Pyomo or CVXOPT.

Solution 2:

Following is a python module that includes a function LPmi(.) to solve mixed integer linear programs. It employs the Branch and Bound algorithm on top of scipy.optimize.linprog(.). It is this responder's creation; anyone is free to use or modify it. It also includes an example in the form of a test(.) function.

import numpy as np
from scipy.optimize import linprog
import copy

class LP:
    minimise = True
    c = None
    A_ub = None
    b_ub = None
    A_eq = None
    b_eq = None
    bounds = None
    method = ""
    fun = 0.
    x = None
    success = False
    
    def __init__(self,c,minimise=True,Aub=None, bub=None, Aeq=None, beq=None,
                 bounds=None,method="revised simplex"):
        self.minimise = minimise
        if self.minimise:
            self.c = c
        else:
            self.c = -1 * c
        self.A_ub = Aub
        self.b_ub = bub
        self.A_eq = Aeq
        self.b_eq = beq
        self.bounds = bounds
        self.method = method
        
    def cal(self): 
        res = linprog(self.c,A_ub=self.A_ub, b_ub=self.b_ub,
                      A_eq=self.A_eq, b_eq=self.b_eq,bounds=self.bounds,
                      method=self.method)                      
        if res["success"] == False:
            return res["message"]
        else:
            self.success = True
        
            if self.minimise:
                self.fun = res["fun"]
            else:
                self.fun = -1 * res["fun"]

            self.x = res["x"]
            
    def get_res(self):
        return (self.x, self.fun, self.success)

class LP_Data:
    
    minimise = True
    c = None
    A_ub = None
    b_ub = None
    A_eq = None
    b_eq = None
    bounds = None
    method = ""
    
    def __init__(self,c,minimise=True,Aub=None, bub=None, Aeq=None, beq=None,
                 bounds=None,method="revised simplex"):
        self.minimise = minimise
        if self.minimise:
            self.c = c
        else:
            self.c = -1 * c
        self.A_ub = Aub
        self.b_ub = bub
        self.A_eq = Aeq
        self.b_eq = beq
        self.bounds = bounds
        self.method = method

    def set_bounds(self,bounds_list):
        self.bounds = bounds_list

class LPd:
    
    data = []
    fun = 0.
    x = []
    success = False
    result = None
    
    def __init__(self,lp_data):
        self.data = lp_data
        self.clip_bound = np.array(list(zip(*self.data.bounds)))

    def cal(self):
        result = None
        res = linprog(self.data.c,A_ub=self.data.A_ub, b_ub=self.data.b_ub,
                      A_eq=self.data.A_eq, b_eq=self.data.b_eq,
                      bounds=self.data.bounds,
                      method=self.data.method)                      
        if res["success"] == False:
            self.result = ([], np.NaN, False, None)
        else:
            self.success = True
            if self.data.minimise:
                self.fun = res["fun"]
            else:
                self.fun = -1 * res["fun"]
            self.x = res["x"].clip(self.clip_bound[0], self.clip_bound[1])
            self.result = (self.x, self.fun, self.success, self.data)
            
    def get_res(self):
        return self.result

def level_iterator(level0, int_index):                                  

    level1 = []
    for k in range(len(level0)):
        for i in int_index:
            if level0[k][0][i-1] != np.floor(level0[k][0][i-1]):
                cur_bounds = level0[k][3].bounds[i-1]
                lp = LPd(copy.deepcopy(level0[k][3]))
                lp.data.bounds[i-1] = (cur_bounds[0],
                       max(cur_bounds[0], floor(level0[k][0][i-1])))
                lp.cal()
                output = lp.get_res()
                level1.append(output)
                lp = LPd(copy.deepcopy(level0[k][3]))
                lp.data.bounds[i-1] = (min(np.ceil(level0[k][0][i-1]),
                      cur_bounds[1]), cur_bounds[1])
                lp.cal()
                output = lp.get_res()
                level1.append(output)
                break
    return level1

def intsol(solution,int_index):
    is_int = True 
    for i in int_index: 
        if solution[0][i-1] != np.floor(solution[0][i-1]):
            is_int = False
            break
    return is_int
 
def feasiblelevl(level, solution, int_index):
    newlevel = []
    solution = solution
    for k in range(len(level)):
        lp = level[k]
        if len(lp[0]) > 0:
            if lp[2] == False:
                pass
            elif intsol(lp,int_index) and lp[1] >= solution[1]:
                solution = lp
            elif lp[1] > solution[1]:
                newlevel.append(lp)
    return (newlevel, solution)

def LPmi(data, int_index):
    level0 = []
    lp  = LPd(data)
    lp.cal()
    level0.append(lp.get_res())
    solution = [None,-np.inf,None,None]
    level0 = (level0, solution)
    level0 = feasiblelevl(level0[0], solution, int_index)
    if len(level0[0]) == 0:
        return level0[1][0:3]
    else:
        for k in range(10):
            level1 = level_iterator(level0[0],int_index)
            level1 = feasiblelevl(level1, solution, int_index)
            if len(level1[0]) == 0:
                break
            level0 = level1
        return level1[1][0:3]
    
def test():
    c = np.array([3.,7])
    minimise = False
    A_ub = np.array([[7.,8],[1,3]])
    b_ub = np.array([56,12])
    data = LP_Data(c,minimise,A_ub,b_ub,bounds=[(0,None),(0,None)])
    int_index = [2] #or [1,2] or [] or [1]#
    out = LPmi(data,int_index)
    print(out)
    
if __name__ == "__main__":
    np.set_printoptions(precision=3,suppress=True)
    test()