在这方面有一个很好的讨论 here .
这里有两个问题。(i) 如何找到一个LP的多个最优解，以及(ii)为什么你要这样做。
首先回答(ii)--通常情况下，你可能想要求所有的最优解决方案，因为有一些次要的或较低的重要性的目标要在它们之间进行选择（例如，想要从当前设置中最小化变化，或以某种方式最小化风险）。如果是这种情况，我的个人建议是找到一种方法，将低阶偏好纳入你的目标函数（例如，通过添加一个低权重的术语）。
在回答(i)时，我发现用图形来观察LP是有帮助的--你所举的例子在这方面很有效（两个变量意味着你可以绘制它--见 关于沃尔弗拉姆的情节 ).你的每个约束条件都在这个不等式图上画了一条线，解决方案只能选择在这条线的一边。
你的目标函数就像在这个可行的区域内有一个恒定的梯度，你试图找到最高的地方。你可以通过将你的目标函数设置为一个特定的值并画出该线来画出目标函数的轮廓。如果你这样做，你会发现你的目标函数轮廓线与顶部的约束线（你的第一个约束线）是平行的。
You can see this directly from the equation: 6x1 + 9x2 <= 100 divides down to 2x1 + 3x2 <= 100/3, and your objective divides down to have the same gradient. What this means is you can move along that top constraint from one corner to the other without changing the value of your objective function.
有无限多的最优解可以解决这个方程。
2x1+3x2==100/3，在x1==0，和x1==20/3之间。你已经确定了两个角的解决方案。
如果你想找到所有同样最优的节点--对于大型问题来说，可能有大量这样的节点--那么下面的代码给出了一个基本的实现方法。 here .当你第一次运行它时，它会给你一个角落的解决方案 - 然后你需要把这个节点（一组变量和懈怠为零）添加到A中，并迭代直到目标下降。你可以把这个放在一个循环中。请注意，在目前的实现中，这只适用于变量下限为0，上限为无界的问题。
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])