线性规划的单纯形法及python实现

1. 线性规划的解

一个线性规划（LP）问题如果有最优解，必定有一个最优解出现在顶点上。

从几何上讲，线性规划的约束集合为凸面体，一个超平面（目标函数的等值面）靠近约束集合时，最先接触到的必然是凸面体的一个顶点或一条边，因此一定有一个最优解出现在顶点上。

2. 线性规划的标准型

线性规划的标准型为：
$$\begin{array}{rl}\min & {\mathbf{c}}^T\mathbf{x}\\ s.t.& \mathbf{Ax}=\mathbf{b}\\ & \mathbf{x}\ge 0\end{array}$$
可以通过以下两个方式转化成标准型：

1. 将所有不等式约束转化为等式约束：通过增加非负的松弛变量，
2. 将所有变量转化为非负约束：负约束变量通过引入非负变量，无约束变量使用两个非负变量的差值表示。

对于线性方程组$\mathbf{Ax}=\mathbf{b}$，$\mathbf{A}\in {\mathbb{R}}^{m\times n}$，通过令$n-m$个变量为0（非基变量），解出的变量称为基变量，构成一组基本解。

3. 初级单纯形法

单纯形法就是在顶点上搜索LP的最优解。求解步骤如下：

1. 初始化：化成标准型，写出单纯形表，并找一个初始基本解。
   例如（见《运筹学（原书第2版）》（——罗纳德L.拉丁）5.3节）：考虑如下标准型：
   $$\begin{array}{rl}\max & 12x_1+9x_2\\ s.t.& x_1+x_3=1000\\ & x_2+x_4=1500\\ & x_1+x_2+x_5=1750\\ & 4x_1+2x_2+x_6=4800\\ & x_1,x_2,x_3,x_4,x_5,x_6\ge 0\end{array}$$
   先找到一个顶点，即初始基本解。令$x_1,x_2=0$为非基变量，可以求解方程组得到初始基本解${\mathbf{x}}^{(0)}=(0,0,1000,1500,1750,4800)$。系数表如下：

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\mathbf{c}$	12	9	0	0	0	0	$\mathbf{b}$
$\mathbf{A}$	1	0	1	0	0	0	1000
	0	1	0	1	0	0	1500
	1	1	0	0	1	0	1750
	4	2	0	0	0	1	4800
	$N$	$N$	$B$	$B$	$B$	$B$
${\mathbf{x}}^{(0)}$	0	0	1000	1500	1750	4800

2. 找到可行方向$\Delta \mathbf{x}$

首先因为相邻顶点的紧约束组只相差一个，并且标准型只有等式约束，而等式约束始终是紧约束，所以这里能放松的只有非基变量的非负约束。我们可以有$(\Delta x_1,\Delta x_2)=(1,0),(0,1)$两种选择。

接着根据线性规划的可行方向的性质可知，可行方向满足$\mathbf{A}\Delta \mathbf{x}=0$，将上面两种情况分别代入，可以求解一个线性方程算出基变量的方向。

然后由$\Delta f={\mathbf{c}}^T\Delta \mathbf{x}$可以算出每个方向对函数值的优化情况。

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\mathbf{c}$	12	9	0	0	0	0	$\mathbf{b}$
$\mathbf{A}$	1	0	1	0	0	0	1000
	0	1	0	1	0	0	1500
	1	1	0	0	1	0	1750
	4	2	0	0	0	1	4800
	$N$	$N$	$B$	$B$	$B$	$B$
${\mathbf{x}}^{(0)}$	0	0	1000	1500	1750	4800	${\mathbf{c}}^T\mathbf{x}=0$
$x_1$的$\Delta \mathbf{x}$	1	0	-1	0	-1	-4	$\Delta f=12>0$
$x_2$的$\Delta \mathbf{x}$	0	1	0	-1	-1	-2	$\Delta f=9>0$

3. 判断最优性
   根据每个方向的$\Delta f$可知函数值还可以继续优化，不能优化则迭代停止。
4. 确定步长
   根据可行方向有负值，可知约束集在此处有界。如果找到一个可以无限移动的方向，说明约束集无界，停止迭代。
   这里取$x_1$的$\Delta \mathbf{x}$作为迭代方向。现在让${\mathbf{x}}^{(0)}$点沿着可行方向移动，直到有一个基变量最先降为0，因此步长为$\lambda =1000$。
5. 新顶点和基
   新顶点为$$\begin{array}{rl}{\mathbf{x}}^{(1)}& ={\mathbf{x}}^{(0)}+\lambda \Delta \mathbf{x}\\ & =(0,0,1000,1500,1750,4800)+1000\times (1,0,-1,0,-1,-4)\\ & =(1000,0,0,1500,750,800)\end{array}$$
   此时，$x_1$入基成为新的基变量，$x_3$出基成为非基变量，新基为 { x 1 , x 4 , x 5 , x 6 } \{ x_1, x_4, x_5, x_6\} {x1​,x4​,x5​,x6​}。回到第二步，继续寻找可行方向，判断是否能继续优化。迭代过程如下：

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\mathbf{c}$	12	9	0	0	0	0	$\mathbf{b}$
$\mathbf{A}$	1	0	1	0	0	0	1000
	0	1	0	1	0	0	1500
	1	1	0	0	1	0	1750
	4	2	0	0	0	1	4800
$t=0$	$N$	$N$	$B$	$B$	$B$	$B$
${\mathbf{x}}^{(0)}$	0	0	1000	1500	1750	4800	${\mathbf{c}}^T\mathbf{x}=0$
$x_1$的$\Delta \mathbf{x}$	1	0	-1	0	-1	-4	$\Delta f=12>0,\lambda =1000$
$x_2$的$\Delta \mathbf{x}$	0	1	0	-1	-1	-2	$\Delta f=9>0$
$t=1$	$B$	$N$	$N$	$B$	$B$	$B$
${\mathbf{x}}^{(1)}$	1000	0	0	1500	750	800	${\mathbf{c}}^T\mathbf{x}=12000$
$x_2$的$\Delta \mathbf{x}$	0	1	0	-1	-1	-2	$\Delta f=9>0,\lambda =400$
$x_3$的$\Delta \mathbf{x}$	-1	0	1	0	0	4	$\Delta f=-12<0$
$t=2$	$B$	$B$	$N$	$B$	$B$	$N$
${\mathbf{x}}^{(2)}$	1000	400	0	1100	350	0	${\mathbf{c}}^T\mathbf{x}=15600$
$x_3$的$\Delta \mathbf{x}$	-1	2	1	-2	-1	0	$\Delta f=6>0,\lambda =350$
$x_6$的$\Delta \mathbf{x}$	0	-0.5	0	0.5	0.5	1	$\Delta f=-4.5<0$
$t=3$	$B$	$B$	$B$	$B$	$N$	$N$
${\mathbf{x}}^{(3)}$	650	1100	350	400	0	0	${\mathbf{c}}^T\mathbf{x}=17700$
$x_5$的$\Delta \mathbf{x}$	1	-2	-1	2	1	0	$\Delta f=-6<0$
$x_6$的$\Delta \mathbf{x}$	-0.5	0.5	0.5	-0.5	0	1	$\Delta f=-1.5<0$

因此最优解为${\mathbf{x}}^{\ast }=(650,1100,350,400,0,0)$。

单纯形法只在两种情况下停止迭代：一是证明该标准型可行域无界，二是找到全局最优解。

4. 修正单纯形法

由于在单纯形法中每次迭代都要求解方程组，为了计算高效，使用矩阵运算进行求解。将基变量的系数列向量组合成矩阵$\mathbf{B}$，那么基向量可以通过${\mathbf{B}}^{-1}\mathbf{b}$得到，单纯形方向通过$-{\mathbf{B}}^{-1}{a}^{(j)}$得到，其中$a^{(j)}$是非基变量$x_j$的系数列向量。

换基之后，矩阵$\mathbf{B}$的某一列被新基的系数列向量所替代，为了快速计算逆矩阵，新基的逆矩阵可以由旧的逆矩阵变换得到：${\mathbf{B}}_{t+1}^{-1}={\mathbf{EB}}_t^{-1}$，其中枢轴（pivot）矩阵
$$\mathbf{E}=\left[\begin{array}{cccccccc}1& 0& \cdots & 0& -\frac{\Delta x_{1st}}{\Delta x_{leave}}& 0& \cdots & 0\\ 0& 1& \cdots & 0& -\frac{\Delta x_{2nd}}{\Delta x_{leave}}& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& & ⋮\\ 0& 0& \cdots & 1& -\frac{\Delta x_{jth}}{\Delta x_{leave}}& 0& \cdots & 0\\ 0& 0& \cdots & 0& -\frac{1}{\Delta x_{leave}}& 0& \cdots & 0\\ 0& 0& \cdots & 0& ⋮& 0& \ddots & 0\\ 0& 0& \cdots & 0& -\frac{\Delta x_{mth}}{\Delta x_{leave}}& 0& \cdots & 1\end{array}\right]$$
是一个类似于单位阵的矩阵，区别在于出基的那一列被一个列向量替换，分母$\Delta x_{leave}$是单纯形方向中出基变量对应的值，分子$\Delta x_{jth}$是单纯形方向中基变量$j$对应的元素，从上到下依次排列，主对角线上的元素为$-1/\Delta x_{leave}$。

有了逆矩阵之后，可以直接算出目标的变化量：
$$\begin{array}{rl}\Delta f& ={\mathbf{c}}^T\Delta \mathbf{x}=c_j+(c_{1st},\dots ,c_{mth})\Delta {\mathbf{x}}^B\\ & =c_j-(c_{1st},\dots ,c_{mth}){\mathbf{B}}^{-1}{a}^{(j)}\\ & \triangleq c_j-\mathbf{v}{a}^{(j)}\end{array}$$
，其中$\Delta {\mathbf{x}}^B$表示现有基的方向，倒数第二个等式使用了前文公式，$\mathbf{v}\triangleq (c_{1st},\dots ,c_{mth}){\mathbf{B}}^{-1}$称为定价向量。因此，无需计算出单纯形方向就可以得到目标值的变化情况。

修正单纯形法的步骤：

1. 初始化。找到一组初始可行基，使用其系数矩阵$\mathbf{B}$求解方程组${\mathbf{Bx}}^B=\mathbf{b}$得到一组初始可行解${\mathbf{x}}^{(0)}$，其中非基变量都是0。
   例如：同样是第3节的例子，选择$(x_3,x_4,x_5,x_6)$为初始基，则系数矩阵如下：

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\mathbf{c}$	12	9	0	0	0	0	$\mathbf{b}$
$\mathbf{A}$	1	0	1	0	0	0	1000
	0	1	0	1	0	0	1500
	1	1	0	0	1	0	1750
	4	2	0	0	0	1	4800
	$N$	$N$	$B$	$B$	$B$	$B$
${\mathbf{x}}^{(0)}$	0	0	1000	1500	1750	4800

$${\mathbf{B}}_{t=0}^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$
3. 定价。使用系数矩阵求出定价向量$\mathbf{v}=({\mathbf{c}}^B{)}^T{\mathbf{B}}^{-1}$，其中${\mathbf{c}}^B$是基变量在目标函数中的系数向量。然后分别令每一个非基变量$x_j=1$，使用公式(3)可以求出每个非基变量的目标函数变化量。
$$\begin{array}{rl}& \mathbf{v}=\left[\begin{array}{cccc}0& 0& 0& 0\end{array}\right],\\ & \\ & \text{when }{x}_1=1,\Delta f=c_1=12,\\ & \text{when }{x}_2=1,\Delta f=c_1=9\\ & \end{array}$$
5. 判断最优解。如果没有一个方向可以继续优化目标函数，则停止迭代，否则选择一个非基变量入基。（注意此时无法判断哪个变量出基，需要在下一步计算单纯形方向后才能决定）

根据$\Delta f$可知$x_1,x_2=1$都能优化目标值，取$x_1=1$。

6. 单纯形方向。$\Delta {\mathbf{x}}^B=-{\mathbf{B}}^{-1}{a}^{(j)}$，$a^{(j)}$是上一步选择入基的变量的系数向量。
   $$\Delta {\mathbf{x}}^B=-{\mathbf{B}}^{-1}{a}^{(1)}=-\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}1\\ 0\\ 1\\ 4\end{array}\right]=\left[\begin{array}{c}-1\\ -0\\ -1\\ -4\end{array}\right]$$
7. 步长。如果单纯形方向所有元素都非负，则停止计算，约束无界。否则选择最先降为0的基变量出基，相应的移动距离为步长。

$\lambda =1000$

8. 新基和新顶点。新顶点为${\mathbf{x}}^{(t+1)}={\mathbf{x}}^{(t)}+\lambda \Delta \mathbf{x}$，构造枢轴矩阵$\mathbf{E}$，并更新系数矩阵的逆矩阵为${\mathbf{EB}}^{-1}$。重复步骤2。

根据单纯形方向可知$x_3$出基，因此
$$\mathbf{E}=\left[\begin{array}{cccc}-\frac{1}{-1}& 0& 0& 0\\ 0& 1& 0& 0\\ -\frac{-1}{-1}& 0& 1& 0\\ -\frac{-4}{-1}& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ -1& 0& 1& 0\\ -4& 0& 0& 1\end{array}\right],{\mathbf{B}}_{t=1}^{-1}=\mathbf{E}{\mathbf{B}}_{t=0}^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ -1& 0& 1& 0\\ -4& 0& 0& 1\end{array}\right]$$

迭代过程如下：

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\mathbf{c}$	12	9	0	0	0	0	$\mathbf{b}$
$\mathbf{A}$	1	0	1	0	0	0	1000
	0	1	0	1	0	0	1500
	1	1	0	0	1	0	1750
	4	2	0	0	0	1	4800
$t=0$	$N$	$N$	$1st$	$2nd$	$3rd$	$4th$
${\mathbf{x}}^{(0)}$	0	0	1000	1500	1750	4800

$${\mathbf{B}}_{t=0}^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],{\mathbf{v}}^T=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\Delta f$	12	9	0	0	0	0
$x_1$的$\Delta \mathbf{x}$	1	0	-1	0	-1	-4	$\lambda =1000$
$t=1$	$1st$	$N$	$N$	$2nd$	$3rd$	$4th$
${\mathbf{x}}^{(1)}$	1000	0	0	1500	750	800	${\mathbf{c}}^T\mathbf{x}=12000$

$$\mathbf{E}=\left[\begin{array}{cccc}-\frac{1}{-1}& 0& 0& 0\\ 0& 1& 0& 0\\ -\frac{-1}{-1}& 0& 1& 0\\ -\frac{-4}{-1}& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ -1& 0& 1& 0\\ -4& 0& 0& 1\end{array}\right],{\mathbf{B}}_{t=1}^{-1}=\mathbf{E}{\mathbf{B}}_{t=0}^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ -1& 0& 1& 0\\ -4& 0& 0& 1\end{array}\right],{\mathbf{v}}^T=\left[\begin{array}{c}12\\ 0\\ 0\\ 0\end{array}\right]$$

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\Delta f$	0	9	-12	0	0	0
$x_2$的$\Delta \mathbf{x}$	0	1	0	-1	-1	-2	$\lambda =400$
$t=2$	$1st$	$4th$	$N$	$2nd$	$3rd$	$N$
${\mathbf{x}}^{(2)}$	1000	400	0	1100	350	0	${\mathbf{c}}^T\mathbf{x}=15600$

$$\begin{array}{rl}& \mathbf{E}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& -0.5\\ 0& 0& 1& -0.5\\ 0& 0& 0& 0.5\end{array}\right],{\mathbf{B}}_{t=2}^{-1}=\mathbf{E}{\mathbf{B}}_{t=1}^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 2& 1& 0& -0.5\\ 1& 0& 1& -0.5\\ -2& 0& 0& 0.5\end{array}\right],\\ & {\mathbf{v}}^T=\left[\begin{array}{cccc}12& 0& 0& 9\end{array}\right]{\mathbf{B}}_{t=2}^{-1}=\left[\begin{array}{c}-6\\ 0\\ 0\\ 4.5\end{array}\right]\end{array}$$

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$
$\Delta f$	0	0	6	0	0	-4.5
$x_3$的$\Delta \mathbf{x}$	-1	2	1	-2	-1	0	$\lambda =350$
$t=3$	$1st$	$4th$	$3rd$	$2nd$	$N$	$N$
${\mathbf{x}}^{(3)}$	650	1100	350	400	0	0	${\mathbf{c}}^T\mathbf{x}=17700$

$$\begin{array}{rl}& \mathbf{E}=\left[\begin{array}{cccc}1& 0& -1& 0\\ 0& 1& -2& 0\\ 0& 0& -1& 0\\ 0& 0& 2& 1\end{array}\right],{\mathbf{B}}_{t=3}^{-1}=\mathbf{E}{\mathbf{B}}_{t=2}^{-1}=\left[\begin{array}{cccc}0& 0& -1& 0.5\\ 0& 1& -2& 0.5\\ 1& 0& 1& -0.5\\ 0& 0& 2& -0.5\end{array}\right],\\ & {\mathbf{v}}^T=\left[\begin{array}{cccc}12& 0& 0& 9\end{array}\right]{\mathbf{B}}_{t=3}^{-1}=\left[\begin{array}{c}0\\ 0\\ 6\\ 1.5\end{array}\right]\end{array}$$

$\Delta f$	0	0	0	0	-6	-1.5	最优

5. 修正单纯形的代码实现

import numpy as np
import pandas as pd
import copy

class Simplex:
    def __init__(self,A, b, c):
        columns=["b"]
        for i in range(A.shape[1]):
            columns.append("x"+str(i+1))
        index = ["c"] + columns[A.shape[1]-A.shape[0]+1:]

        b = np.asmatrix(b)
        matrix = np.concatenate((b.T, A), axis=1)
        c = np.asmatrix(np.append(0, c))
        matrix = np.concatenate((c, matrix), axis=0)
        matrix = pd.DataFrame(matrix,index=index, columns=columns)
        self.matrix = copy.deepcopy(matrix)
        # 基变量
        self.basic_varible = self.matrix.iloc[1:,0].index.tolist()
        self.basic_varible_len = len(self.basic_varible)
        # 非基变量
        self.nonbasic_varible = list(set(self.matrix.columns.tolist()[1:])-set(self.basic_varible))
        self.nonbasic_varible_len = len(self.nonbasic_varible)
        # 对目标函数的优化量
        self.reduce_cost = self.matrix.iloc[0][self.nonbasic_varible] 
        print(matrix)

    def solve(self):
        iteration = 1
        B = self.matrix.loc[self.basic_varible,self.basic_varible]
        c = copy.deepcopy(self.matrix.iloc[0,1:])
        while self.reduce_cost.min() < -1e-5:
            print(f"第{iteration}次迭代")
            #出基，进基
            #选择对函数值改进量最大的非基变量xj
            print("目标改进量：\n", self.reduce_cost)
            basic_in = self.reduce_cost.idxmin()
            #非基变量xj的系数aj
            aj = self.matrix.iloc[1:][basic_in]
            deltax = -B.dot(aj) # 单纯形方向
            print("单纯形方向：\n", deltax)

            if deltax.min() > 0:
                print("无界")
                return 

            xB = self.matrix.iloc[1:,0] # 基变量的值
            if xB.min() < 0:
                print("无解")
                return
            
            # 选出基变量
            step = (-xB/deltax).replace([np.inf, -np.inf], np.nan)
            step = step[step>0]
            # print("每个基变量的步长：\n", step)
            lambda_ = step.min()
            basic_out = step.idxmin()
            print(f"进基: {basic_in}",f"出基: {basic_out}")

            self.basic_varible[self.basic_varible.index(basic_out)] = basic_in
            self.nonbasic_varible[self.nonbasic_varible.index(basic_in)] = basic_out

            # 枢轴矩阵
            E = pd.DataFrame(np.eye(self.basic_varible_len), index=B.index,columns=B.columns)
            E.loc[:, basic_out] = -deltax/deltax.loc[basic_out]
            E.loc[basic_out,basic_out] = -1/deltax.loc[basic_out]
            B = E.dot(B)
            print("枢轴矩阵:\n", E)
            
            B.index = self.basic_varible
            B.columns = self.basic_varible
            print("系数矩阵:\n", B)

            cB = c[self.basic_varible]
            v = cB.T.dot(B)
            index1 = self.matrix.iloc[0:,0].index.tolist()
            index1[1:] = self.basic_varible
            self.matrix.index = index1

            self.matrix.iloc[0,0] = 0
            for idx in self.basic_varible:
                if idx == basic_in:
                    self.matrix.loc[idx, self.matrix.columns[0]] = lambda_
                else:
                    self.matrix.loc[idx,self.matrix.columns[0]] += lambda_*deltax.loc[idx]
                self.matrix.iloc[0,0] += self.matrix.loc[idx,:][0] * c[idx]

            for idx in self.matrix.columns[1:]:
                self.matrix.loc[self.matrix.index[0],idx] = c.loc[idx] - v.T.dot(self.matrix.loc[:, idx][1:])
            self.reduce_cost = self.matrix.iloc[0,1:]
            iteration += 1
            print(self.matrix)

        print(f"*********运行完毕，一共迭代{iteration-1}次**************")
        print("******最优解**********")
        for x in self.matrix.columns[1:]:
            if x in self.basic_varible:
                print(f"{x}:",self.matrix.loc[x,"b"],end="\t")
            else:
                print(f"{x}:",0,end="\t")
        print("\nobj:",-self.matrix.iloc[0][0])



A=np.asarray([
    [1, 0, 1, 0, 0 ,0],
    [0, 1, 0, 1, 0, 0],
    [ 1, 1, 0, 0, 1, 0],
    [4, 2, 0, 0, 0, 1]])
b = np.asarray([1000,1500,1750,4800])
c = np.asarray([-12,-9,0,0,0,0])
simplex = Simplex(A,b,c)

simplex.solve()