1. 补充问题
上一节中的代码在运行时还有很多细节没有处理，这里补充两个比较重要的情况：
- 存在等式约束
如果有等式约束，那么就没法通过添加松弛变量直接给出初始可行解，需要用大M法或者两阶段法求初始可行解。计算机求解一般使用二阶段法，即首先给等式约束条件添加人工变量，使得问题有一个初始可行解。注意人工变量和松弛变量/剩余变量不一样。松弛变量/剩余变量是用于将不等式变为等式（标准形），而人工变量则是在等式约束中额外添加一个变量，这个变量最后一定要等于0才行，否则就是没有初始可行解，不能进入第二阶段。
在第一阶段，每个不等式约束对应一个松弛变量，每个等式约束对应一个人工变量。另外，b<0的不等式约束也需要添加人工变量，因为对应的松弛变量等于b不能作为基变量。 - 没有有限最优解（目标函数无界）
对应的情况是：目标函数中有某个系数大于0，对应变量在约束条件中的系数全都≤0，问题结束，没有最优解。
2. python算法实现
首先假设问题除了不等式约束，还有等式约束：
求解问题为：min c*x
s.t.A_ub * x ≤ b_ub
A_eq * x = b_eq
定义maxiter为最大迭代次数，tol为求解精度，OptimizeResult为格式输出函数，_solve_simplex为上一节的求解函数。完整的求解步骤如下：
def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, maxiter=1000, tol=1.0E-12,):
cc = np.asarray(c)
f0 = 0
n = len(c)
mub = len(b_ub)
meq = len(b_eq)
m = mub+meq
n_slack = mub
n_artificial = meq + np.count_nonzero(bub < 0)
Aub_rows, Aub_cols = A_ub.shape
Aeq_rows, Aeq_cols = A_eq.shape
slcount = 0
avcount = 0
status = 0
messages = {0: "Optimization terminated successfully.",
1: "Iteration limit reached.",
2: "Optimization failed. Unable to find a feasible"
" starting point.",
3: "Optimization failed. The problem appears to be unbounded.",
4: "Optimization failed. Singular matrix encountered."}
# 建立第一阶段单纯型表。
T = np.zeros([m+2, n+n_slack+n_artificial+1])
T[-2, :n] = cc
T[-2, -1] = f0
b = T[:-2, -1]
if meq > 0:
T[:meq, :n] = Aeq
b[:meq] = beq
if mub > 0:
T[meq:meq+mub, :n] = Aub
b[meq:meq+mub] = bub
np.fill_diagonal(T[meq:m, n:n+n_slack], 1)
#第一阶段的基向量basis，记录下标。
basis = np.zeros(m, dtype=int)
r_artificial = np.zeros(n_artificial, dtype=int)
for i in range(m):
if i < meq or b[i] < 0:
basis[i] = n+n_slack+avcount
r_artificial[avcount] = i
avcount += 1
if b[i] < 0:
b[i] *= -1
T[i, :-1] *= -1
T[i, basis[i]] = 1
T[-1, basis[i]] = 1
else:
basis[i] = n+slcount
slcount += 1
for r in r_artificial:
T[-1, :] = T[-1, :] - T[r, :]
#第一阶段求解
nit1, status = _solve_simplex(T, n, basis)
# 如果第一阶段的目标函数为0，则删除人工变量，进入第二阶段
if abs(T[-1, -1]) <= sol:
T = T[:-1, :]
T = np.delete(T, np.s_[n+n_slack:n+n_slack+n_artificial], 1)
else:
status = 2
if status != 0:
message = messages[status]
if disp:
print(message)
return OptimizeResult(x=np.nan, fun=-T[-1, -1], nit=nit1, status=status, message=message, success=False)
nit2, status = _solve_simplex(T, n, basis, maxiter=maxiter-nit1, phase=2, callback=callback, tol=tol, nit0=nit1)
solution = np.zeros(n+n_slack+n_artificial)
solution[basis[:m]] = T[:m, -1]
x = solution[:n]
slack = solution[n:n+n_slack]
obj = -T[-1, -1]
return OptimizeResult(x=x, fun=obj, nit=int(nit2), status=status, slack=slack, message=messages[status], success=(status == 0))
下面是例子：
import numpy as np
c=np.array([2,3,-5])
A_ub=np.array([[-2,5,-1],[1,3,1]])
B_ub=np.array([-10,12])
A_eq=np.array([[1,1,1]])
B_eq=np.array([7])
res=linprog(-c,A_ub,B_ub,A_eq,B_eq)
print(res)
输出为
fun: -14.571428571428571
message: 'Optimization terminated successfully.'
nit: 2
slack: array([3.85714286, 0.57142857, 6.42857143, 7. , 0. ])
status: 0
success: True
x: array([6.42857143, 0.57142857, 0. ])
