Добро пожаловать в SO!
Здесь есть хорошее обсуждение этого .
Здесь есть два вопроса: (i) Как найти несколько оптимальных решений для LP, и (ii) Почему вы хотели бы это сделать.
Сначала ответьте на (ii) - обычно вы можете запросить все оптимальные решения, потому что между ними стоит выбрать какую-то второстепенную или более низкую цель (для пример, желающий минимизировать изменение текущего параметра или каким-либо образом минимизировать риск). Если это так, моя личная рекомендация состояла бы в том, чтобы найти способ включить это предпочтение более низкого порядка в вашу целевую функцию (например, добавив член с низким весом).
Ответ (i) Я считаю это помогает посмотреть на LP графически - и приведенный вами пример отлично подходит для этого (две переменные означают, что вы можете построить его - см. график на вольфраме ). Каждое из ваших ограничений помещает линии на этот график неравенства, и решения могут быть выбраны только на одной стороне этой линии. 
Ваша целевая функция похожа на постоянный градиент в этой возможной области, и вы пытаетесь найти самую высокую точку. Вы можете нарисовать контуры своей целевой функции, установив для нее определенное значение c и нарисовав эту линию. Если вы это сделаете, вы обнаружите, что контуры вашей целевой функции параллельны верхней линии ограничения (ваше первое ограничение).
Вы можете увидеть это непосредственно из уравнения: 6x1 + 9x2 <= 100 делит до 2x1 + 3x2 <= 100/3, и ваша цель делится, чтобы иметь тот же градиент. Это означает, что вы можете перемещаться по этому верхнему ограничению из одного угла в другой, не меняя значения вашей целевой функции. </p>
Существует бесконечно много оптимальных решений, которые решают уравнение:
2x1 + 3x2 == 100/3, между x1 == 0 и x1 == 20/3. Вы уже определили решения по двум углам.
Если вы хотите найти все узлы, которые одинаково оптимальны - для больших проблем их может быть большое количество - тогда приведенный ниже код дает основную c реализация описанного здесь метода здесь . Когда вы запустите его в первый раз, он даст вам одно из угловых решений - затем вам нужно добавить этот узел (набор переменных и нулевых резервов) в A и повторять до тех пор, пока цель не ухудшится. Вы можете поместить это в al oop.
import pulp as pulp
# Accounting:
# n structural varuables (n = 2)
# m constraints (m = 2)
# => No. of basics = 2 (no. of constraints)
# => No. of non-basics = 2 (no. of variables)
nb = 2
M = 100 # large M value - upper bound for x1, x2 * the slacks
model = pulp.LpProblem('get all basis', pulp.LpMaximize)
# Variables
x = pulp.LpVariable.dicts('x', range(2), lowBound=0, upBound=None, cat='Continuous')
# Non-negative Slack Variables - one for each constraint
s = pulp.LpVariable.dicts('s', range(2), lowBound=0, upBound=None, cat='Continuous')
# Basis variables (binary)
# one for each variable & one for each constraint (& so slack)
B_x = pulp.LpVariable.dicts('b_x', range(len(x)), cat='Binary')
B_s = pulp.LpVariable.dicts('b_s', range(len(s)), cat='Binary')
# Objective
model += 2000*x[0] + 3000*x[1]
# Constraints - with explicit slacks
model += 6*x[0] + 9*x[1] + s[0] == 100
model += 2*x[0] + x[1] + s[1] == 20
# No. of basics is correct:
model += pulp.lpSum(B_x) + pulp.lpSum(B_s) == nb
# Enforce basic and non-basic behaviour
for i in range(len(x)):
model += x[i] <= M*B_x[i]
for i in range(len(x)):
model += s[i] <= M*B_s[i]
# Cuts - already discovered solutions
A = []
# A = [[1, 1, 0, 0]]
# A = [[1, 1, 0, 0], [0, 1, 0, 1]]
for a in A:
model += (B_x[0]*a[0] + B_x[1]*a[1] +
B_s[0]*a[2] + B_s[1]*a[3]) <= nb - 1
model.solve()
print('Status:', pulp.LpStatus[model.status])
print('Objective:', pulp.value(model.objective))
for v in model.variables():
print (v.name, "=", v.varValue)