运筹系列1：线性规划单纯形法python代码

1. 线性规划问题

线性规划标准模型如下：

$\min cx$
s.t. $Ax=b$
$x\ge 0$

首先来看最简单的例子：A是一个方阵，比如：
$\min 2x_1+5x_2$
s.t. $3x_1+3x_2=6$
$4x_1+x_2=5$
$x_1,x_2\ge 0$

只有一个确定解，最优值为7。

约束条件不够的时候，表现为A变扁了。这时候的技巧就是，将A竖着切一刀，分离出一个方阵来：

基变量：令$x=(x_B,x_N)$，其中$x_B$满秩，称$x_B$称为基变量，$x_N$称为非基变量。

将原问题用基变量重新表达为

$\min c_B{x}_B+c_N{x}_N$
s.t. $B{x}_B+N{x}_N=b$
$x_B\ge 0,x_N\ge 0$

极点：令$x$有$m$维（秩$m_B\le m$），则需要$m$个等式才能求解，它们的交点称为极点（有时候是极线）。

单纯形法将搜索问题限制在了极点/极线范围内。那么如何找到它们呢？如有$x_N=0$，则$x$就是极点了。那么此时有：

主问题可行条件：$B^{-1}b\ge 0$

主问题极点有很多，哪一个才是最优解呢？除了遍历搜索，还有一种方法：看对偶问题是否可行。上面问题的对偶问题为：

$\max by$
s.t. $y{B}^T\le c_B^T$
$y{N}^T\le c_N^T$

B转置后仍然是满秩的方阵，因此等号是可以取到的，极点满足$y{B}^T=c_B^T$，那么此时有：

对偶问题可行条件：$c_B^T{B}^{-1}N\le c_N^T$

如果我们基向量拆分出来同时满足上面两个可行条件，那么恭喜找到最优解了：

主问题可行：非基变量=0时，基变量系数$B^{-1}b\ge 0$；
对偶问题可行：基变量系数=0时，非基变量系数$c_N^T-c_B^T{B}^{-1}N\ge 0$.
定理：线性规划的极点同时满足主问题和对偶问题的可行条件时，这个极点是最优解。

这两个条件非常重要！重要！重要！

2. 求解步骤

定理：线性规划问题的最优解一定是极点或极线，并且满足比相邻极点都要优。

单纯形法的步骤是从一个初始极点出发，不断找到更优的相邻极点，直到找到最优的极点（或极线）。

消去$x_B$得到问题的字典表达，即：
$\min c_B^T{B}^{-1}b+(c_N^T-c_B^T{B}^{-1}N)x_N$
s.t. $B^{-1}b-B^{-1}N{x}_N\ge 0$
$x_N\ge 0$

称$(c_N^T-c_B^T{B}^{-1}N{)}_i$为第$i$个非基变量的残差(reduced cost)，如果残差小于0，那么说明这个非基变量还没有满足最优条件。直观上来看，这个非基变量取稍稍大于0的数就可以继续优化目标函数了。我们选择一个基变量和这个非基变量对换，就可以找到更优的相邻极点。
另一种理解方法：残差对应的是对偶问题可行条件，小于0表示对偶问题还不可行，需要继续探索。

接下来的问题是具体选择哪个基变量和哪个非基变量进行对换。我们用启发式原则，每次将负数残差$(c_N^T-c_B^T{B}^{-1}N{)}_i$最小（绝对值最大）的非基变量$x_i$替换为基变量，同时将$(\frac{B^{-1}b}{(B^{-1}N{)}_i}{)}_j$最小值对应的基变量$x_j$替换为非基变量。这个进基/出基的过程称为pivoting。
另一种表达方式是：$\min z=cx$，s.t.$Ax=b$的pivoting是每次找出c中最小数对应的非基变量$x_i$，再找出$b_i/A_{ij}$最小的基变量$x_j$进行对换。

如果是max问题，令$c^{\prime }=-c$即可转化为min问题，相对应的，每次pivoting是找出$c^{\prime }$最大值对应的非基变量。

3. python算法实现

这里假设原问题都是小于等于约束，这样添加松弛变量之后，问题一定有初始可行解；同时假设问题存在有限最优解。特殊情况将在下一节进行处理。代码为：

import numpy as np   
def solve(d):
    (bn,cn) = d.shape
    s = list(range(cn-bn,cn-1)) #基变量列表
    while max(d[-1][:-1]) > 0:
        jnum = np.argmax(d[-1][:-1]) #转入下标
        out = d[:-1,-1]/np.maximum(0,d[:-1,jnum])
        np.place(out,np.isnan(out),np.inf) # np.nan替换为np.inf
        inum = np.argmin(out)  #转出下标
        s[inum] = jnum #更新基变量
        d[inum]/=d[inum][jnum]
        for i in range(bn):
            if i != inum:
                d[i] -= d[i][jnum] * d[inum]
    return d,s
            
def printSol(d,s):
    for i in range(d.shape[1] - 1):
        print("x%d=%.2f" % (i,d[s.index(i)][-1] if i in s else 0))
    print("objective is %.2f"%(-d[-1][-1]))

求解问题：max $z=x_0+14\ast x_1+6\ast x_2$
s.t.
$x_0+x_1+x_2\le 4$
$x_0\le 2$
$x_2\le 3$
$3\ast x_1+x_2\le 6$

添加变量将问题变为标准形，保存到data.txt中：

1 1 1 1 0 0 0 4
1 0 0 0 1 0 0 2
0 0 1 0 0 1 0 3
0 3 1 0 0 0 1 6
1 14 6 0 0 0 0 0

调用代码：

d = np.loadtxt("data.txt", dtype=np.float)
d,s = solve(d)
printSol(d,s)

或者直接从字符串解析：

dstr = '''1 1 1 1 0 0 0 4
1 0 0 0 1 0 0 2
0 0 1 0 0 1 0 3
0 3 1 0 0 0 1 6
1 14 6 0 0 0 0 0'''
d = np.array([[eval(dij) for dij in dj.split(' ')] for dj in dstr.splitlines()]).astype(np.float)
d,s = solve(d)
printSol(d,s)

运行后输出结果为：

x0=0.00
x1=1.00
x2=3.00
x3=0.00
x4=2.00
x5=0.00
x6=0.00
objective is 32.00

4. 代码跟踪

我们尝试跟踪下评论区的一个例子：

初始化：

第一次pivoting：

第二次pivoting：

第三次pivoting：

此时检验数行（也就是最后一行）都是负数了，因此得到结果：
‘’’
x0=22.53
x1=23.20
x2=7.33
x3=0.00
x4=0.00
x5=0.00
objective is 135.27
‘’’

5. 写后感

将simplex用代码写出来，才觉得以前纠结那么久的问题原来那么简单。两三行代码能说清楚的事，何必写一堆看得人眼花缭乱的数学公式呢。
另外，线性规划还有一些很基础的理论要掌握好：

1. 极点和极方向的理论，这个是单纯型法的理论基础。可以参考这里
2. 对偶理论，这个在以后经常会用到。