单纯形法原理

单纯形法

参考连接：单纯形法

单纯形法是针对求解线性规划问题的一个算法，这个名称里的 “单纯形” 是代数拓扑里的一个概念，可以简单将“单纯形”理解为一个凸集，标准的线性规划问题（线性规划标准型）可以表示为：

\[max\,(or\,min)\quad{f(x)=cx}\\ \quad\quad\quad{s.t.}\quad{Ax=b}\\ \quad\quad\quad\quad\quad{x}\geq{0} \]

单纯形法最早由 George Dantzig于1947年提出，单纯形法对于求解线性规划问题是具有跨时代意义的，其实不仅仅是针对线性规划，非线性规划问题在求解的过程中也大量依赖单纯形法。

一、凸集及其性质的介绍

1.1 凸集的概念

单纯形可以理解为一个凸集，介绍单纯形法之前有必要先来了解一下凸集概念及其性质，凸集上求最优解是凸优化的一个分支，凸集对于简化问题是有重要意义的。

凸集可以用图形表示为一个没有洞的联通区域，如下图所示：

凸集的概念：

• （1）如果集合\(C\)中任意两个点\(X_1,X_2\)，其连线上所有的点也都是集合\(C\)中的点，称\(C\)为凸集

• （2）\(X_1,X_2\)的连线可以表示为：\(aX_1+(1-a)X_2\quad(0<a<1)\)

• （3）凸集定义用数学解析式可表示为：对\(\forall{X_1,X_2}\in{C}\)，有：\(aX_1+(1-a)X_2\in{C}\quad(0<a<1)\)

关于（2），以一个一维坐标例子简单理解一下。如图，\(aX_1+(1-a)X_2=a(X_1-X_2)+X_2\quad(0<a<1)\)，其中\(a=0\)和\(a=1\)分别代表\(X_2、X_1\)；\(a\in(0,1)\)代表点\(X_1\)和点\(X_2\)连线中任意一点

根据定义下面这个图形所定义的区域不是凸集：

1.2 凸集的顶点

顶点：

凸集\(C\)中满足下述条件的点\(X\)称为顶点：

• 如果\(C\)中不存在任何两个不同的点\(X_1,X_2\)，使\(X\)成为这两个点连线上的一个点。

• 或者这样表述：对\(\forall{X_1,X_2}\in{C}\)，不存在\(X=aX_1+(1-a)X_2\quad(0<a<1)\)，则称\(X\)是凸集\(C\)的顶点

如图，因为\(a≠0或1\)，所以\(X\)不能成为\(C\)中某条连线上的一个点。

1.3 几个基本定理

定理1    若线性规划问题存在可行解，则问题的可行域是凸集

• 证明见《运筹学教程（第五版）》第20页

定理2    线性规划问题的基可行解\(X\)对应线性规划问题可行域（凸集）的顶点

定理3    若线性规划问题有最优解，一定存在一个基可行解是最优解

二、单纯形法原理

2.1 线性规划问题的解的概念

• 可行解    满足所有约束条件的解，称为线性规划问题的可行解。全部可行解的集合称为可行域。

• 最优解    使目标函数达到最大值的可行解称为最优解。

• 基    设\(A\)为约束方程组的\(m×n(m>n)\)阶系数矩阵，其秩为\(m\)，\(B\)是\(A\)的一个\(m×m\)阶满秩子矩阵，称\(B\)是线性规划问题的一个基。

  • \(B\)中每一个列向量\(P_j(j=1,\dots,m)\)称为基向量

  • 与基向量\(P_j\)对应的变量\(x_j\)称为基变量

  • 线性规划中除基变量以外的变量称为非基变量

• 基解    由\(m\)个约束方程可解出\(m\)个基变量的唯一解\(X_B\)，将这个解加上非基变量取\(0\)的值产生的解\(X\)，称为线性规划问题的基解。

  • 基解中变量取非\(0\)值的个数不大于方程数\(m\)，故基解的总数不超过\(C^m_n\)个

• 基可行解    满足变量非负约束条件（如\(x\geq0\)）的基解称为基可行解

  • 单纯型法应用于标准化形式的线性规划问题，所以在线性规划问题中决策变量\(x\)均大于等于\(0\)

• 可行基    对应于基可行解的基称为可行基

2.2 单纯形法的实现

根据定理2和定理3，最优解一定是一个基可行解，且为凸集的一个顶点。

单纯形法的核心思想可以归纳为： 找到每一个基本可行解，代入目标函数后计算函数值取其最大或最小值即可。单纯形法上从代数角度是寻找约束条件的每一个基本可行解，从几何意义上来说是遍历凸集的每一个顶点，根据算法的特性有时也称为转轴法。

当求解最小值问题时，如果借助梯度的概念，在得到一个顶点后，切换下一个顶点方向是函数变小的方向无疑会节省很多时间。

2.2.1 单纯形法的迭代

约束条件系数矩阵\(A\)中，将作为基的列集合记作\(B\)（或者说\(B\)为线性规划问题的一个基）；非基的列集合记作\(N\)，则有：\(A=(N,B)\)。其中，\(N、B\)均为\(A\)的一个分块矩阵。

同样地，把\(x=\left[\begin{matrix}x_{1}\\x_{2}\\\vdots\\x_n\end{matrix}\right]\)分量中基变量的的集合记作\(x_B\)，非基变量集合记为\(x_N\)

约束条件可写为：\((N,B)(x_N,x_B)=b\)，得到等式\(Nx_N+Bx_b=b\)，矩阵\(B\)是基矩阵且可逆，得到：

\[Bx_B=b-Nx_N\\ B^{-1}Bx_B=B^{-1}b-B^{-1}Nx_N\\ x_B=B^{-1}b-B^{-1}Nx_N\tag{1} \]

当\(x_N=0\)时，\(x_B=B^{-1}b\)，若\(B^{-1}b\geq0\)，则称\(x_B\)为初始基可行解

同样地，将目标函数中的系数向量\(c\)也分为两部分，基变量对应的系数向量记为\(c_B\)，非基变量对应的系数记为\(c_N\)

由此，可得到目标函数的另一种表达形式：

\[f(x)=c_Bx_B+c_Nx_N\tag{2} \]

把(1)代入(2)中，得到目标函数另一种表达形式：

\[f(x)=c_B(B^{-1}b-B^{-1}Nx_N)+c_Nx_N=c_BB^{-1}b-(c_BB^{-1}N-c_N)x_N\tag{3} \]

若此时存在初始基可行解，即\(x_N=0，x_B=B^{-1}b\)，代入(3)后可得到：

\[f(x)=c_BB^{-1}b\tag{4} \]

对(3)进行求导，可得：

\[\frac{\partial{f}}{\partial{x}}=-(c_BB^{-1}N-c_N)\tag{5} \]

若\(\frac{\partial{f}}{\partial{x}}<0\)，即\(c_BB^{-1}N-c_N>0\)，说明此时目标函数是单调递减的，也就是在取初始基可行解后目标函数还有变小的空间。此时可以选择矩阵中某个非基变量（非基集合中的某一列）来替代现有基变量（替代现有基），即选择入基。

如何选择入基呢？函数的微分形式有：

\[f(x+\Delta{x})=f(x)+\frac{\partial{f}}{\partial{x}}\Delta{x}\tag{6} \]

• 补充说明：\(\frac{\partial{f}}{\partial{x}}=\frac{\Delta{f(x)}}{\Delta{x}}\)，所以：\(\Delta{f(x)}=\frac{\partial{f}}{\partial{x}}\Delta{x}\)。又因为：\(\Delta{f(x)}=f(x+\Delta{x})-f(x)\)，所以\(f(x+\Delta{x})=f(x)+\Delta{f(x)}=f(x)+\frac{\partial{f}}{\partial{x}}\Delta{x}\)

由(5)可得，\(c_BB^{-1}N-c_N\)越大则\(\frac{\partial{f}}{\partial{x}}\)越小，目标函数变小的幅度也就越大。将此作为选择入基的标准，遍历每个非基列\(p_j\)（\(p_j\)为\(N\)中的一个列向量），由于每次都只选择一个列作为入基，\(x_N\)其他分量为\(0\)

所以，\(c_BB^{-1}N-c_N\)简化为\(c_BB^{-1}p_j-c_j\)，引入每一个待选列检验量\(\sigma_j=c_BB^{-1}p_j-c_j\)，如果检验量大于\(0\)，就选择这其中检验量最大的列作为入基就能保证目标函数在最大程度上变小，如果待选列检验量都小于等于\(0\)，则代表此时已经得到最小值。

---

入基已经准备就绪，那么如何选择出基呢？即如何选择被替代的基列？在选择入基是其实还有一个隐含的问题，入基被选择后，新的基变量值是多少？这些可以从公式(1)中得到答案：

\[x_B=B^{-1}b-B^{-1}Nx_N=B^{-1}b-B^{-1}p_jx_j\tag{7} \]

\(x_j\)原来是\(x_N\)中分量，初始时为0，将这个非基变量相关列选择作为新的基，\(x_j\)也将获得新的值。

这里设\(\widetilde{b}=B^{-1}b,\,y_j=B^{-1}p_j\)。很明显，\(\widetilde{b}\)是初始基可行解，可行解\(x_B\)可写成矩阵形式：

\[x_B=\left[\begin{matrix}\widetilde{b_1}\\\widetilde{b_2}\\\vdots\\\widetilde{b_m}\end{matrix}\right] -\left[\begin{matrix}y_{j1}\\y_{j2}\\\vdots\\y_{jm}\end{matrix}\right]x_j =\left[\begin{matrix}x_{B1}\\x_{B2}\\\vdots\\x_{Bm}\end{matrix}\right] -\left[\begin{matrix}y_{j1}\\y_{j2}\\\vdots\\y_{jm}\end{matrix}\right]x_j =\left[\begin{matrix}x_{B1}-y_{j1}x_j\\x_{B2}-y_{j2}x_j\\\vdots\\x_{Bm}-y_{jm}x_j\end{matrix}\right]\tag{8} \]

由于约束条件\(x_B\geq0\)的限制，所以\(x_j\)的范围也有限制，即\(x_{Bi}-y_{ji}x_j\geq0\)，也就是\(x_j\leq\frac{x_{Bi}}{y_{ji}}\)

故此，当\(x_j\)取\(min\{\frac{x_{Bi}}{y_{ji}}\}\quad(y_{ji}>0)\)时，即可保证新生成的\(x_B\geq0\)

得到\(x_j\)的值后代入（8），使得原\(x_B\)中有一个基变量变为\(0\)，我们知道\(x_B\)中每一个分量都是与基的列一一对应的，如果其中一个分量变量变为0，则代表这个基被替换，分量变为0的基为出基。

2.2.3 单纯形法代码

\[min\quad{-4x_1-x_2}\\ s.t.\quad-x_1+2x_2\leq4\\ \quad\quad2x_1+3x_2\leq12\\ \quad\quad{x_1}-x_2\leq3\\ \quad\quad{x_1},x_2\geq0 \]

import numpy as np
BASEINDEX=2
#列变换实现单纯形法
def simplex_ColTranslate(c ,A,b,flagstable):B=A[:,BASEINDEX:]Binv=np.linalg.inv(B)xb=np.dot(Binv  ,b)cb=c[:,BASEINDEX:]f=np.dot(cb,xb )Inbaseindex=-1OutBaseIndex = -1tempvalue=0for i in range(BASEINDEX):v=np.squeeze(np.dot(np.dot(cb,Binv),A[:,i])-c[:,i],0)if tempvalue<=v:#找出入基tempvalue=vInbaseindex=iif tempvalue==0:#找到最优解return  xb,f,c,flagstableelse:#入基出基替换，同时交换系数cy=np.dot(Binv,A[:,Inbaseindex])minivalue=Nonefor i in range(y.shape[0]):if y[i]>0:x=xb[i]/y[ i]#找出最小值，保证xb>=0if minivalue is None or minivalue>x :minivalue=xOutBaseIndex=i+BASEINDEXfor i in range(A.shape[1]):if i==OutBaseIndex:for j in range(A.shape[0]):tmp=A[j,i].copy()A[j,i]=A[j, Inbaseindex]A[j, Inbaseindex]=tmptempc= np.squeeze( c[:,OutBaseIndex],0).copy()c[:,OutBaseIndex]=c[:,Inbaseindex]c[:,Inbaseindex]=tempcflagstable[OutBaseIndex-BASEINDEX]=Inbaseindexreturn simplex_ColTranslate(c, A, b,flagstable) def useColTranslate():BASEINDEX = 2c = np.array([[-4, -1, 0, 0, 0]],dtype=np.float)A = np.array([[-1, 2, 1, 0, 0], [2, 3, 0, 1, 0], [1, -1, 0, 0, 1]],dtype=np.float)b = np.array([[4], [12], [3]],dtype=np.float)# 记录下标变动情况flagstable = {}for i in range(c.shape[1] - BASEINDEX):flagstable[i] = BASEINDEX + ix, f, c, flags = simplex_ColTranslate(c, A, b, flagstable)print('最优解:%.2f' % (f[0, 0]))for j in range(x.shape[0]):print('x%d=%.1f' % (flags[j] + 1, x[j, 0])) if __name__ == '__main__':useColTranslate()