1. 補(bǔ)充問題
上一節(jié)中的代碼在運(yùn)行時(shí)還有很多細(xì)節(jié)沒有處理亚铁，這里補(bǔ)充兩個(gè)比較重要的情況：
- 存在等式約束
如果有等式約束，那么就沒法通過添加松弛變量直接給出初始可行解，需要用大M法或者兩階段法求初始可行解。計(jì)算機(jī)求解一般使用二階段法偎行，即首先給等式約束條件添加人工變量，使得問題有一個(gè)初始可行解。注意人工變量和松弛變量/剩余變量不一樣蛤袒。松弛變量/剩余變量是用于將不等式變?yōu)榈仁剑?biāo)準(zhǔn)形）熄云，而人工變量則是在等式約束中額外添加一個(gè)變量，這個(gè)變量最后一定要等于0才行汗盘，否則就是沒有初始可行解皱碘，不能進(jìn)入第二階段。
在第一階段隐孽，每個(gè)不等式約束對應(yīng)一個(gè)松弛變量癌椿，每個(gè)等式約束對應(yīng)一個(gè)人工變量。另外菱阵，b<0的不等式約束也需要添加人工變量踢俄，因?yàn)閷?yīng)的松弛變量等于b不能作為基變量。 - 沒有有限最優(yōu)解（目標(biāo)函數(shù)無界）
對應(yīng)的情況是：目標(biāo)函數(shù)中有某個(gè)系數(shù)大于0晴及，對應(yīng)變量在約束條件中的系數(shù)全都≤0都办，問題結(jié)束，沒有最優(yōu)解虑稼。
2. python算法實(shí)現(xiàn)
首先假設(shè)問題除了不等式約束琳钉，還有等式約束：
求解問題為：min c*x
s.t.A_ub * x ≤ b_ub
A_eq * x = b_eq
定義maxiter為最大迭代次數(shù)，tol為求解精度蛛倦，OptimizeResult為格式輸出函數(shù)歌懒，_solve_simplex為上一節(jié)的求解函數(shù)。完整的求解步驟如下：
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及皂，記錄下標(biāo)。
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)
# 如果第一階段的目標(biāo)函數(shù)為0且改，則刪除人工變量验烧，進(jìn)入第二階段
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. ])
