基础不牢，地动山摇，致敬很多年前夏铱铱的一句话。
相信大家看了这么多期精彩的推文，也大概能推测出我的工作。言归正传，在列生成一文中提起过单纯形法以及最优性检验以木材切割为例讲解列生成法Column Generation(附Python代码）。那会用的是求解器Gurobi处理限制主问题（一个线性规划问题），单纯形法是求解线性规划问题的基础。在上运筹学课程时，一般在单纯形表上迭代计算线性规划问题，今天就写一个代码来代替人工计算，暂时讨论有最优解或无界的情况，并且只考虑松弛变量（不考虑剩余变量（很多情况下有人把它和松弛变量混为一谈，其实它两还是有区别的），人工变量），退化啥的都暂时不考虑。
存在最优解时的终止条件
到这里已经知道，当所有非基变量的Reduce cost都大于等于0的时候，最优解找到了。
进基与出基
进基：选择Reduce cost最小并且小于0, 所对应的非基变量xj，作为这一次迭代的进基。
出基：取进基xj所在列的系数向量aj，令 sita = min{bi/aji | aji >0} = bk/ajk,则xk为出基。
这里顺带说一下，什么情况下会出现问题无界，存在非基变量的Reduce cost小于0，但是它所在列的系数向量小于等于0，则无界，证明略。
在确定好进基和出基之后，进行一个矩阵的行初等变换，单纯形表。假设单纯形表是一个m*n的矩阵S。(1)xk所在行每一个元素除以S[xk,xj],（2）完成(1)后，S[xk,xj] = 1，将xj所在列的系数除了S[xk,xj]外的其他元素通过行变换变为0,（3）在基变量中，用xj代替xk。重复（1）（2）（3），直至最优或无界。
案例分析
一个最小化线性规划问题
min -x1-x2
x1 + (9/14)x2 <= 51/14
-2x1 + x2 <= 1/3
x1,x2 >= 0
将其标准化
min -x1-x2 + 0*x3 + 0*x4
x1 + (9/14)x2 +x3 == 51/14
-2x1 + x2 +x4== 1/3
x1,x2 ,x3,x4>= 0
则初始单纯形表如下，第一行对应着变量的Reduce cost, S[obj,b],这个格子是目标函数的相反数，基变量x3 = 51/14 ,x4=1/3 ,非基变量x1,x2,Reduce cost 分别为 -1，-1，那就选下标最小的非基变量x1作为进基好了，a1 = [1.0,-2.0],b=[3.64,0.33], 很简单就知道x3是出基，看上面的出基选取规则。
第一次迭代，还有Reduce cost 小于0，继续迭代。
第二次迭代，好了，Reduce cost 都大于0了，最优解找到了。
最优解和目标函数值，嗯，obj忘记取相反数了，尴尬！
完整代码
import numpy as npimport pandas as pdimport copyclass Simplex:def __init__(self,matrix):self.matrix = copy.deepcopy(matrix)self.basic_varible = self.matrix.iloc[1:,0].index.tolist()self.nonbasic_varible = list(set(self.matrix.columns.tolist()[1:])-set(self.basic_varible))self.reduce_cost = self.matrix.iloc[0,1:]print("初始单纯形表")print(matrix)def solve(self):iteration = 1while self.reduce_cost.min() < 0:print(f"第{iteration}次迭代")#出基，进基basic_out,basic_in = self.select_basic()print(f"进基: {basic_in}",f"出基: {basic_out}")#更新单纯形表#(1)将basic_out所在行除以self.matrix.loc[basic_out,basic_in]self.matrix.loc[basic_out,:] = self.matrix.loc[basic_out,:]/self.matrix.loc[basic_out,basic_in]#(2)矩阵的初等行变换，basic_in所在列，除了self.matrix.loc[basic_out,basic_in]是1，其他都是0for idx in self.matrix.index:if idx != basic_out:self.matrix.loc[idx,:] = self.matrix.loc[idx,:] -self.matrix.loc[idx,basic_in]*self.matrix.loc[basic_out,:]#(3)在单纯形表的基变量中，用进基替换出基index1 = self.matrix.index.tolist()i = index1.index(basic_out)index1[i] = basic_inself.matrix.index = index1#更新基变量与非基变量self.basic_varible = self.matrix.iloc[1:,0].index.tolist()self.nonbasic_varible = list(set(self.matrix.columns.tolist()[1:])-set(self.basic_varible))print(self.matrix)self.reduce_cost = self.matrix.iloc[0,1:]iteration += 1print(f"*********运行完毕，一共迭代{iteration-1}次**************")print("******最优解**********")for x in self.basic_varible:print(f"{x}:",self.matrix.loc[x,"b"],end="\t")print("\nobj:",-self.matrix.iloc[0][0])def select_basic(self):#选择reduce cost最小的非基变量xjbasic_in = self.reduce_cost.idxmin()#非基变量xj的系数P = self.matrix.iloc[1:][basic_in]if P.max() <= 0:print("无界")returnb = self.matrix.iloc[1:,0]if b.min() < 0:print("无解")returntmp = P.index[0]i = 0while P[tmp] <= 0:m += 1tmp = P.index[m]for idx in P.index:if P[idx] > 0:if b[idx]/P[idx] < b[tmp]/P[tmp]:tmp = idxbasic_out = tmpreturn basic_out,basic_in
命名都挺朴实无华的，对照着前面的解说，应该能看懂。其实是备案给我自己以后复习的。
运行
matrix = pd.DataFrame(np.array([[0,-1,-1,0,0],[51/14,1,9/14,1,0],[1/3,-2,1,0,1]]),index=["obj","x3","x4"],columns=["b","x1","x2","x3","x4"])simplex = Simplex(matrix)simplex.solve()
运行结果
初始单纯形表b x1 x2 x3 x4obj 0.000000 -1.0 -1.000000 0.0 0.0x3 3.642857 1.0 0.642857 1.0 0.0x4 0.333333 -2.0 1.000000 0.0 1.0第1次迭代进基: x1 出基: x3b x1 x2 x3 x4obj 3.642857 0.0 -0.357143 1.0 0.0x1 3.642857 1.0 0.642857 1.0 0.0x4 7.619048 0.0 2.285714 2.0 1.0第2次迭代进基: x2 出基: x4b x1 x2 x3 x4obj 4.833333 0.0 0.0 1.3125 0.15625x1 1.500000 1.0 0.0 0.4375 -0.28125x2 3.333333 0.0 1.0 0.8750 0.43750*********运行完毕，一共迭代2次********************最优解**********x1: 1.4999999999999996 x2: 3.3333333333333335obj: -4.833333333333333
OK，一个完整的解说已经完成。如果不是出于学习目的，还是乖乖趁着还有学生身份，申请求解器的学术许可吧。后期可能是，对偶单纯形法和割平面，记得关注点赞。
