How to get all combinations other than objective value using python and Linear programming

Question

I obtained the max profit of 6442143.99530000 and the quota for each product by using pulp. However, I would like to output the next suboptimal solutions other than the optimal solution in rank order. Is it possible? In other words, I want to get all combinations assigned to each product.

I need the help of several experts here. Thanks in advance for your help.

[temp.csv]

product,profit,min_alloc,max_alloc
Prod_A,19251.1503,5,50
Prod_B,11029.62628,5,50
Prod_C,49455.97051,5,50
Prod_D,2270.916437,5,50
Prod_E,20727.92984,5,50
Prod_F,41979.71183,5,50
Prod_G,10298.78459,5,50
Prod_H,9912.822331,5,50
Prod_I,4579.744941,5,50
Prod_J,25079.03564,5,50
Prod_K,3754.784861,5,50

[np_test.py]

import pandas as pd
import numpy as np
from pulp import *

pd.options.display.float_format = '{:.5f}'.format
total_capa = 200

df = pd.read_csv('temp.csv')

profit_sum = sum(df['profit'])

LP = LpProblem(name = "LP", sense = LpMaximize)

X1 = LpVariable(name='Prod_A', lowBound=5, upBound=50, cat='Integer')
X2 = LpVariable(name='Prod_B', lowBound=5, upBound=50, cat='Integer')
X3 = LpVariable(name='Prod_C', lowBound=5, upBound=50, cat='Integer')
X4 = LpVariable(name='Prod_D', lowBound=5, upBound=50, cat='Integer')
X5 = LpVariable(name='Prod_E', lowBound=5, upBound=50, cat='Integer')
X6 = LpVariable(name='Prod_F', lowBound=5, upBound=50, cat='Integer')
X7 = LpVariable(name='Prod_G', lowBound=5, upBound=50, cat='Integer')
X8 = LpVariable(name='Prod_H', lowBound=5, upBound=50, cat='Integer')
X9 = LpVariable(name='Prod_I', lowBound=5, upBound=50, cat='Integer')
X10 = LpVariable(name='Prod_J', lowBound=5, upBound=50, cat='Integer')
X11 = LpVariable(name='Prod_K', lowBound=5, upBound=50, cat='Integer')

np_product = np.array([X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11])


for profit in df['profit']:
    np_profit = np.array(df['profit'])
    

# OBJECTIVE function
LP.objective = sum(np_profit * np_product)
print('Objective: ', LP.objective, '\n')
   
# CONSTRAINTS
constraints = [
    X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 + X11 == total_capa,
]

print('Constraints: ')
for i in constraints:
    print(i)

for i, c in enumerate(constraints):
    constraint_name = f"const_{i}"
    LP.constraints[constraint_name] = c

# SOLVE model
res = LP.solve()

# results
print('\nResults: ')
for v in LP.variables():
    print(f"{v} Product: {v.varValue:5.1f} EA")

print('\n', 'The objective funtion is ', value(LP.objective))

Result - Optimal solution found

Objective value:                6442143.99530000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01


Results:
Prod_A Product:   5.0 EA
Prod_B Product:   5.0 EA
Prod_C Product:  50.0 EA
Prod_D Product:   5.0 EA
Prod_E Product:  15.0 EA
Prod_F Product:  50.0 EA
Prod_G Product:   5.0 EA
Prod_H Product:   5.0 EA
Prod_I Product:   5.0 EA
Prod_J Product:  50.0 EA
Prod_K Product:   5.0 EA

 The objective funtion is  6442143.9953

Answer 1

AFAIK, there are two ways to achieve this:

1. Use a commercial solver like Gurobi and its solution pool to collect the N best solutions.

2. Solve your model, add an integer cut to exclude the found solution and solve the model again. Repeat until you have enough solutions. One way to model such an integer cut is explained here.

After cleaning up your code a bit, here's a working example to find the next 5 best solutions:

import pandas as pd
import numpy as np
import pulp


def print_solution(x):
    print('\nResults: ')
    for x in x.values():
        print(f"{x.name}: {x.varValue:5.1f} EA")
    print('\n', 'The objective funtion is ', pulp.value(prob.objective))


df = pd.read_csv('temp.csv')
profit = df["profit"].values
total_capa = 200

# the problem
prob = pulp.LpProblem(name="LP", sense=pulp.LpMaximize)

# Add variables
x = {}
for i, name in enumerate(df["product"].values):
    x[i] = pulp.LpVariable(name=name, lowBound=5, upBound=50, cat="Integer")

# OBJECTIVE function
prob.objective = sum(profit[i] * x[i] for i in range(len(x)))
print('Objective: ', prob.objective, '\n')

# CONSTRAINTS
prob.addConstraint(pulp.lpSum(x[i] for i in range(len(x))) == total_capa)

# SOLVE model
res = prob.solve()
print_solution(x)

# Add integer cut to exclude the optimal solution
bigM = 45
num_sols = 5
z = []
b = []
for k in range(num_sols):
    # helper variables
    z.append([pulp.LpVariable(name=f"z[{k},{i},]", cat="Integer") for i in range(len(x))])
    b.append([pulp.LpVariable(name=f"b[{k},{i},]", cat="Binary") for i in range(len(x))])

    prob.addConstraint(pulp.lpSum(z[k][i] for i in range(len(x))) >= 1)
    for i in range(len(x)):
        prob.addConstraint(z[k][i] <= x[i] - x[i].varValue + bigM * b[k][i])
        prob.addConstraint(z[k][i] <= -x[i] + x[i].varValue + bigM * (1-b[k][i]))

    # solve the model again
    prob.solve()
    print_solution(x)