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

运筹学 专栏收录该内容
73 篇文章 80 订阅

1. 线性规划问题

线性规划标准模型如下:

min ⁡ c x \min cx mincx
s.t. A x = b Ax = b Ax=b
x ≥ 0 x\ge 0 x0

首先来看最简单的例子:A是一个方阵,比如:
min ⁡ 2 x 1 + 5 x 2 \min 2x_1+5x_2 min2x1+5x2
s.t. 3 x 1 + 3 x 2 = 6 3x_1+3x_2=6 3x1+3x2=6
4 x 1 + x 2 = 5 4x_1+x_2=5 4x1+x2=5
x 1 , x 2 ≥ 0 x_1,x_2 \ge 0 x1,x20

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

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

基变量:令 x = ( x B , x N ) x=(x_B,x_N) x=(xB,xN),其中 x B x_B xB满秩,称 x B x_B xB称为基变量 x N x_N xN称为非基变量

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

min ⁡ c B x B + c N x N \min c_Bx_B+c_Nx_N mincBxB+cNxN
s.t. B x B + N x N = b Bx_B+Nx_N=b BxB+NxN=b
x B ≥ 0 , x N ≥ 0 x_B\ge 0,x_N\ge 0 xB0,xN0

极点:令 x x x m m m维(秩 m B ≤ m m_B\le m mBm),则需要 m m m个等式才能求解,它们的交点称为极点(有时候是极线)。

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

主问题可行条件: B − 1 b ≥ 0 B^{-1}b\ge 0 B1b0

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

max ⁡ b y \max by maxby
s.t. y B T ≤ c B T yB^T\le c^T_B yBTcBT
y N T ≤ c N T yN^T\le c^T_N yNTcNT

B转置后仍然是满秩的方阵,因此等号是可以取到的,极点满足 y B T = c B T yB^T= c^T_B yBT=cBT,那么此时有:

对偶问题可行条件: c B T B − 1 N ≤ c N T c^T_BB^{-1}N\le c^T_N cBTB1NcNT

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

主问题可行:非基变量=0时,基变量系数 B − 1 b ≥ 0 B^{-1}b\ge 0 B1b0
对偶问题可行:基变量系数=0时,非基变量系数 c N T − c B T B − 1 N ≥ 0 c^T_N-c^T_BB^{-1}N\ge 0 cNTcBTB1N0.
定理:线性规划的极点同时满足主问题和对偶问题的可行条件时,这个极点是最优解。

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

2. 求解步骤

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

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

消去 x B x_B xB得到问题的字典表达,即:
min ⁡ c B T B − 1 b + ( c N T − c B T B − 1 N ) x N \min c^T_BB^{-1}b+(c^T_N-c^T_BB^{-1}N)x_N mincBTB1b+(cNTcBTB1N)xN
s.t. B − 1 b − B − 1 N x N ≥ 0 B^{-1}b-B^{-1}Nx_N\ge 0 B1bB1NxN0
x N ≥ 0 x_N\ge 0 xN0

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

接下来的问题是具体选择哪个基变量和哪个非基变量进行对换。我们用启发式原则,每次将负数残差 ( c N T − c B T B − 1 N ) i (c^T_N-c^T_BB^{-1}N)_i (cNTcBTB1N)i最小(绝对值最大)的非基变量 x i x_i xi替换为基变量,同时将 ( B − 1 b ( B − 1 N ) i ) j (\frac{B^{-1}b}{(B^{-1}N)_i})_j ((B1N)iB1b)j最小值对应的基变量 x j x_j xj替换为非基变量。这个进基/出基的过程称为pivoting
另一种表达方式是: min ⁡ z = c x \min z = cx minz=cx,s.t. A x = b Ax = b Ax=b的pivoting是每次找出c中最小数对应的非基变量 x i x_i xi,再找出 b i / A i j b_i/A_{ij} bi/Aij最小的基变量 x j x_j xj进行对换。

如果是max问题,令 c ′ = − c c'=-c c=c即可转化为min问题,相对应的,每次pivoting是找出 c ′ c' c最大值对应的非基变量。

3. python算法实现

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

import numpy as np   
def solve():
    while max(d[-1][:-1]) > 0:
        jnum = np.argmax(d[-1][:-1]) #转入下标
        inum = np.argmin(d[:-1,-1]/d[:-1,jnum])  #转出下标
        s[inum] = jnum #更新基变量
        d[inum]/=d[inum][jnum]
        for i in range(bn):
            if i != inum:
                d[i] -= d[i][jnum] * d[inum]
            
def printSol():
    for i in range(cn - 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 ∗ x 1 + 6 ∗ x 2 z = x_0+14*x_1+6*x_2 z=x0+14x1+6x2
s.t.
x 0 + x 1 + x 2 ≤ 4 x_0 + x_1 + x_2 \leq 4 x0+x1+x24
x 0 ≤ 2 x_0 \leq 2 x02
x 2 ≤ 3 x_2 \leq3 x23
3 ∗ x 1 + x 2 ≤ 6 3*x_1 + x_2 \leq 6 3x1+x26

添加变量将问题变为标准形,保存到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)
(bn,cn) = d.shape
s = list(range(cn-bn,cn-1)) #基变量列表
solve()
printSol()

运行后输出结果为:

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. 写后感

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

  1. 极点和极方向的理论,这个是单纯型法的理论基础。可以参考这里
  2. 对偶理论,这个在以后经常会用到。
©️2021 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值