线性规划——单纯形法（原理及代码实现）

线性规划基本模型

由于单纯性法是用于求解线性规划模型，因此我们对一般的问题进行松弛之后得到了线性规划模型的一般形式（及是LP问题的一般形式）如下：
$$\begin{gathered} maxz=c^{T}X \\ s.t.\{\begin{array}{l}AX=b\\ X\ge 0\end{array} \end{gathered}$$

单纯形法的矩阵描述

针对于该模型引入基变量$X_B$和非基变量$X_N$、系数矩阵的基$B$和非基部分$N$以及基变量对应的系数数列$c_B^T$、非基变量对应的系数数列$c_N^T$可以将原问题化为如下形式：
$$\begin{gathered} maxz=c_B^T{X}_B+c_N^T{X}_N \\ s.t.\{\begin{array}{l}B{X}_B+N{X}_n=b\\ X_B,X_N\ge 0\end{array} \end{gathered}$$
将目标函数与约束条件进行联立，消去$X_B$可以得到目标函数的新的表达式如下所示：
$$z=C_B{B}^{-1}b+(C_N-C_B{B}^{-1}N)X_N$$
可以发现进行联立之后计算式所得的目标函数与基变量无关，要确定从非基变量进行的检验数为$C_N-C_B{B}^{-1}N$,若想要目标函数变得更大则检验数需要大于零，从而可以确定进基变量为非基变量中检验数大于零的一个。

因此对于单纯行表可以表示为如下：

	基变量$X_B$	非基变量$X_N$	右侧RHS
系数矩阵	$I$	$B^{-1}N$	$B^{-1}b$
检验数	0	$C_N-C_B{B}^{-1}N$	$-C_B{B}^{-1}b$

任意时刻各个部分的核心是某个已知矩阵的部分左乘一个$B^{-1}$，因此求解的核心在于快速地维护 $B^{-1}$。

以下我们设$P_k$ 是 $x_k$ 对应的原始系数矩阵的那一列。
又存在递推式：
$$B_i^{-1}=E_i{B}_{i-1}^{-1}$$
其中$E_i$是把一个单位矩阵中，第$j$列替换为${\xi }_i$后的结果，其中$j$表示本次新换入的基在$B_i$中对应的第$j$列，${\xi }_i$由本次换入变量在换入前$B_{i-1}^{-1}{N}_{i-1}$中对应的列$(a_1,a_2,....,a_m)$变换得到，设$l$是换出变量对应的行，则
$${\xi }_i=(-\frac{a_1}{a_l},...,\frac{1}{a_l},...,-\frac{a_m}{a_l})$$
于是，
$$B_i^{-1}=(e_1,...,e_{j-1},{\xi }_j,e_{j+1}...,e_m)B_{i-1}^{-1}$$

寻找入基变量

入基变量求解根据检验数
$${\sigma }_i=C_{N_i}-C_{B_i}{B}_i^{-1}{N}_i$$
中找最大值下标即可以找到。

寻找出基变量

出基变量根据$\theta$法则求出
$$\theta =min\left\{\frac{B^{-1}b}{B_i^{-1}{P}_k}|B_i^{-1}{P}_k>0\right\}$$
即可以找到。
当所有检验数都小于零则单纯形法结束，因为目标值已经不可能再大一步了。

python实现单纯形法

使用单纯形法对如下线性规划模型进行求解：
$$\begin{gathered} maxz=70x_1+30x_2 \\ s.t.\{\begin{array}{l}3x_1+9x_2+x_3=540\\ 5x_1+5x_2+x_4=450\\ 9x_1+3x_2+x_5=720\\ x_1,x_2,x_3,x_4,x_5\ge 0\end{array} \end{gathered}$$
具体实现代码如下所示：

# 单纯形法
import pandas as pd
import numpy as np
from pandas import DataFrame


# 定义线性规划求解函数
def lp_solver(matrix: DataFrame):
    '''
    输入线性规划的矩阵，根据单纯形法求解线性规划模型
    max cx
    s.t. ax <= b
    矩阵形式是：
        b   x1  x2  x3  X4  X5
    obj 0   70  30  0   0   0
    x3  540 3   9   1   0   0
    x4  450 5   5   0   1   0
    x5  720 9   3   0   0   1
    第1行是目标函数的系数
    第2-4行是约束方程的系数
    第1列是约束方程的常数项
    obj-b 交叉，即第1行第一列的元素是目标函数的负值
    x3,x4,x5 既是松弛变量，也是初始可行解
    :param matrix
    :return
    '''
    # 检验数是否大于零
    c = matrix.iloc[0, 1:]
    while c.max() > 0:
        # 选择入基变量
        c = matrix.iloc[0, 1:]
        in_x = c.idxmax()
        # 入基变量的系数
        in_x_v = c[in_x]
        # 选择出基变量
        # 选择正的最小比值对应的变量出基 min(b列/入基变量列）
        b = matrix.iloc[1:, 0]
        # 选择入基变量对应的列
        in_x_a = matrix.iloc[1:][in_x]
        # 得到出基变量
        out_x = (b / in_x_a).idxmin()
        # 旋转操作
        print(matrix.loc[out_x, in_x])
        matrix.loc[out_x, :] = matrix.loc[out_x, :] / matrix.loc[out_x, in_x]
        for idx in matrix.index:
            if idx != out_x:
                matrix.loc[idx, :] = matrix.loc[idx, :] - matrix.loc[out_x, :] * matrix.loc[idx, in_x]
        # 索引替换（入基与出基变量名称替换）
        index = matrix.index.tolist()
        i = index.index(out_x)
        index[i] = in_x
        matrix.index = index
    # 打印结果
    print("最终的最优单纯形法是：")
    print(matrix)
    print("目标函数值是：", -matrix.iloc[0, 0])
    print("最优决策变量是：")
    x_count = (matrix.shape[1] - 1) - (matrix.shape[0] - 1)
    x = matrix.iloc[0, 1:].index.tolist()[: x_count]
    for xi in x:
        print(xi, '=', matrix.loc[xi, 'b'])


# 主程序代码
if __name__ == '__main__':
    # 约束方程系数矩阵，包含常数项
    matrix = pd.DataFrame(
        np.array([
            [0, 70, 30, 0, 0, 0],
            [540, 3, 9, 1, 0, 0],
            [450, 5, 5, 0, 1, 0],
            [720, 9, 3, 0, 0, 1], ]),
        index=['obj', 'x3', 'x4', 'x5'],
        columns=['b', 'x1', 'x2', 'x3', 'x4', 'x5']
    )
    # 调用前面定义的函数求解
    lp_solver(matrix)

最终的求解结果如下：

最终的最优单纯形法是：
          b   x1   x2   x3   x4        x5
obj -5700.0  0.0  0.0  0.0 -2.0 -6.666667
x3    180.0  0.0  0.0  1.0 -2.4  1.000000
x2     15.0  0.0  1.0  0.0  0.3 -0.166667
x1     75.0  1.0  0.0  0.0 -0.1  0.166667
目标函数值是： 5700.0
最优决策变量是：
x1 = 75.0
x2 = 15.0