Multiple optimal solutions in PuLP

Question

I have an LP that returns one optimal solution using PulP. I suspect there are many optimal solutions to my LP, but I don't know how to get at them. Even getting the two optimal corner points would do and then I could at least say that some linear combination of the two points are optimal.

Here is my LP:

        my_lp_problem = pulp.LpProblem("My LP Problem", pulp.LpMinimize)
        x1=pulp.LpVariable('x1', lowBound=0.1, upBound=95, cat='Continuous')
        x2=pulp.LpVariable('x2', lowBound=0.1, upBound=95, cat='Continuous')
        x3=pulp.LpVariable('x3', lowBound=0.1, upBound=95, cat='Continuous')
        x4=pulp.LpVariable('x4', lowBound=0.1, upBound=95, cat='Continuous')
        x5=pulp.LpVariable('x5', lowBound=0.1, upBound=95, cat='Continuous')
        x6=pulp.LpVariable('x6', lowBound=0.1, upBound=95, cat='Continuous')
        x7=pulp.LpVariable('x7', lowBound=0.1, upBound=95, cat='Continuous')
        x8=pulp.LpVariable('x8', lowBound=0.1, upBound=95, cat='Continuous')
        x9=pulp.LpVariable('x9', lowBound=0.1, upBound=95, cat='Continuous')
        y1=pulp.LpVariable('y1', lowBound=0.1, upBound=95, cat='Continuous')
        y2=pulp.LpVariable('y2', lowBound=0.1, upBound=95, cat='Continuous')
        y3=pulp.LpVariable('y3', lowBound=0.1, upBound=95, cat='Continuous')

        x_prime=pulp.LpVariable('x_prime', lowBound=0, cat='Continuous')
        matrix_sum0=sum(matrix[0])
        matrix_sum1=sum(matrix[1])
        matrix_sum2=sum(matrix[2])
        matrix_sum3=sum(matrix[3])
        matrix_sum4=sum(matrix[4])
        matrix_sum5=sum(matrix[5])
        matrix_sum6=sum(matrix[6])
        matrix_sum7=sum(matrix[7])
        matrix_sum8=sum(matrix[8])
        matrix_sum9=sum(matrix[9])
        matrix_sum10=sum(matrix[10])
        matrix_sum11=sum(matrix[11])

        property_sum=np.sum(weights_per_property)

        my_lp_problem+=x_prime, 'Z'        

        my_lp_problem+=matrix_sum0*x1+matrix_sum1*x2+matrix_sum2*x3+matrix_sum3*x4+matrix_sum4*x5+matrix_sum5*x6+matrix_sum6*x7+matrix_sum7*x8+matrix_sum8*x9+matrix_sum9*y1+matrix_sum10*y2+matrix_sum11*y3-property_sum<=x_prime

        my_lp_problem+=-(matrix_sum0*x1+matrix_sum1*x2+matrix_sum2*x3+matrix_sum3*x4+matrix_sum4*x5+matrix_sum5*x6+matrix_sum6*x7+matrix_sum7*x8+matrix_sum8*x9+matrix_sum9*y1+matrix_sum10*y2+matrix_sum11*y3-property_sum)<=x_prime

        my_lp_problem+=x1+x2+x3+x4+x5+x6+x7+x8+x9+y1+y2+y3==100
        my_lp_problem+=x1>=x2
        my_lp_problem+=x2>=x3
        my_lp_problem+=x3>=x4
        my_lp_problem+=x4>=x5
        my_lp_problem+=x5>=x6
        my_lp_problem+=x6>=x7
        my_lp_problem+=x7>=x8
        my_lp_problem+=x8>=x9
        my_lp_problem+=x9>=y1
        my_lp_problem+=y1>=y2
        my_lp_problem+=y2>=y3

        my_lp_problem+=matrix[0][0]*x1+matrix[1][0]*x2+matrix[2][0]*x3+matrix[3][0]*x4+matrix[4][0]*x5+matrix[5][0]*x6+matrix[6][0]*x7+matrix[7][0]*x8+matrix[8][0]*x9+matrix[9][0]*y1+matrix[10][0]*y2+matrix[11][0]*y3==weights_per_property[0]
        my_lp_problem+=matrix[0][1]*x1+matrix[1][1]*x2+matrix[2][1]*x3+matrix[3][1]*x4+matrix[4][1]*x5+matrix[5][1]*x6+matrix[6][1]*x7+matrix[7][1]*x8+matrix[8][1]*x9+matrix[9][1]*y1+matrix[10][1]*y2+matrix[11][1]*y3==weights_per_property[1]
        my_lp_problem+=matrix[0][2]*x1+matrix[1][2]*x2+matrix[2][2]*x3+matrix[3][2]*x4+matrix[4][2]*x5+matrix[5][2]*x6+matrix[6][2]*x7+matrix[7][2]*x8+matrix[8][2]*x9+matrix[9][2]*y1+matrix[10][2]*y2+matrix[11][2]*y3==weights_per_property[2]
        my_lp_problem+=matrix[0][3]*x1+matrix[1][3]*x2+matrix[2][3]*x3+matrix[3][3]*x4+matrix[4][3]*x5+matrix[5][3]*x6+matrix[6][3]*x7+matrix[7][3]*x8+matrix[8][3]*x9+matrix[9][3]*y1+matrix[10][3]*y2+matrix[11][3]*y3==weights_per_property[3]
        my_lp_problem+=matrix[0][4]*x1+matrix[1][4]*x2+matrix[2][4]*x3+matrix[3][4]*x4+matrix[4][4]*x5+matrix[5][4]*x6+matrix[6][4]*x7+matrix[7][4]*x8+matrix[8][4]*x9+matrix[9][4]*y1+matrix[10][4]*y2+matrix[11][4]*y3==weights_per_property[4]
        my_lp_problem+=matrix[0][5]*x1+matrix[1][5]*x2+matrix[2][5]*x3+matrix[3][5]*x4+matrix[4][5]*x5+matrix[5][5]*x6+matrix[6][5]*x7+matrix[7][5]*x8+matrix[8][5]*x9+matrix[9][5]*y1+matrix[10][5]*y2+matrix[11][5]*y3==weights_per_property[5]
        my_lp_problem+=matrix[0][6]*x1+matrix[1][6]*x2+matrix[2][6]*x3+matrix[3][6]*x4+matrix[4][6]*x5+matrix[5][6]*x6+matrix[6][6]*x7+matrix[7][6]*x8+matrix[8][6]*x9+matrix[9][6]*y1+matrix[10][6]*y2+matrix[11][6]*y3==weights_per_property[6]

        my_lp_problem.solve()
        out_file.write(str([variable.varValue for variable in my_lp_problem.variables()])[1:-1] + '\n')    
        pulp.LpStatus[my_lp_problem.status]
        for variable in my_lp_problem.variables():
            print("{} = {}".format(variable.name, variable.varValue))

        print('objective value = ' + str(pulp.value(my_lp_problem.objective)))
        print(my_lp_problem.status)

The optimal solution to the objective function is 0 and that is what occurs when PuLP provides me with my x1, ..., y3. But I know there are more x1, ..., y3 out there that will also provide a value of 0 for the objective function. I would like to get some of those, perhaps even just the two corner points that give the optimal solution. Is there a way to do this by solving and looping?

Answer 1

Sorry, this is not easy to do using an LP solver. With some effort, we can encode bases using binary variables and then enumerate the optimal bases (i.e. optimal "corner points"). This means solving a series of MIP models. Here is a link with more information.