LP线性规划求解 之 单纯形 算法

LP线性规划求解 之 单纯形 算法

认识-单纯形

• 核心: 顶点旋转

1. 随机找到一个初始的基本可行解

2. 不断沿着可行域旋转(pivot)

3. 重复2,直到结果不能改进为止

案例-过程

以上篇的case2的松弛型为例.

$min \ y = -x1-x2$

s.t.

$50x1 + 20x2 + a1 = 2000 \\ -1.5x1+x2 + a2 =0 \\ x1-x2+a3=0 \\ x1,x2,a1,a2,a3 >=0\\ 其中a1,a2,a3为松弛变量$

即:

基本变量(松弛): a1, a2, a3

非基本变^量: x1, x2

$min \ y=-x1-x2\\ s.t.$

$a1 = 2000-50 x1-20 x2 \\ a2 = 1.5x1 -x2 \\ a3 = -x1+x2 \\ x1,x2,a1,a2,a3 >=0$

基本可行解即将右边非基本变量设置为0, 计算左边基本变量(松弛)的值

令: x1, x2为0, 即基本可行解为:

$x^t = (0,0,2000,0,0)^t, 目标值y= (-x1-x2) = 0$

旋转(pivot)流程

• 替入变量: 当前目标函数中第一个系数为负的非基本变量(松弛)为替入变量

• 替出变量: 最严格约束里的基本变量为替出变量

即:

• x1 是目标函数中第一个系数为负数(替入变量), 计算x1在各条件下的约束值

• case1: x1 = 40

• case2: x1 = $\infty$

• case3: x1 = 0

得到条件3为最严格约束, 即对应的a3为替出变量.

即由 $a3 = -x1+x2 \ 得 x1 = x2-a3$

对原目标及约束中的x1进行替换即:

$min \ y = -2x2+a3 \\ s.t.$

$x1 = x2-a3 \\ a1 = 2000-50(x2-a3) - 20x2 = 2000 - 70x2 + 50a3 \\ a2 = 1.5(x2 - a3) - x2 = 0.5x2 - 1.5a3$

然后旋转,令右边变量为0, 即:

$x = (0,0,0,2000,0)^t, 目标y=0$

发现目标值y没有变小, 于是再继续转

• 目标中第一为负的是x2, 其在各条件下的约束值为

• case1: $x2 =\infty$

• case2: $x2 = \frac {2000}{70}$

• case3: $x2 = \infty$

得到条件2为最严格约束, 即对应的a1为替出变量.

即由$a1 = 2000-70x2+50a3 \ 得 \ x2 = \frac {2000 + 50a3 -a1}{70}$

对原目标的x2作替换即:

$min \ y = \frac{-4000}{70} - \frac{3}{7} a3 + \frac{1}{35}xa1 \\ s.t.$

$x1 =\frac{2000}{70} - \frac{2}{7}a3- \frac{1}{70}a1 \\ x2 = \frac {2000 + 50a3 -a1}{70} \\ a2 = \frac {1000}{70} - \frac {16}{14}a3 - \frac {1}{140}a1$

然后再旋转, 令右边的变量为0, 即:

$x = (\frac {200}{70}, \frac {2000}{70}, 0, \frac{1000}{70},0), 目标 \ y = - \frac{4000}{70}$

然后再继续旋转...

....

当目标函数没有系数为负, 即说明无法再通过旋转继续减小目标函数, 此时收敛,即停止旋转, 此时的基本解即为最优解.

最后得 $x = (25,37,5,12.5,0,0), 达到最优解目标 \ min \ y = -62.5$

单纯形算法步骤(参考网上)

1. 将LP的目标函数及约束转为标准型

2. 转换为松弛型 (即将不等于转为等于)

3. 找到一个基本可行解, 寻找目标的替入变量, 约束的替出变量, 替换目标函数, 不断旋转

4. 重复步骤3, 直到目标系数中没有系数为负数的项, 收敛, 即结束算法

(附) Python代码实现

• 是搬运大佬的代码, 自己有很多也尚未看明白, 只作个笔记, 当然还是为了copy
• 个人觉得这个大佬的代码,追求简洁但严重丢失了可读性, PEP8呢
• 注释 是后面加的, 当然,原来的版本也基本没啥注释, 很是任性

import numpy as np

class Simplex(object):
    def __init__(self, obj, max_mode=False):
        self.max_mode = max_mode  # 默认是求min，如果是max目标函数要乘-1
        self.mat = np.array([[0] + obj]) * (-1 if max_mode else 1)     

    def add_constraint(self, a, b):
        self.mat = np.vstack([self.mat, [b] + a]) #矩阵加入约束

    def solve(self):
        m, n = self.mat.shape  # 矩阵里有1行目标函数，m - 1行约束，应加入m-1个松弛变量
        temp, B = np.vstack([np.zeros((1, m - 1)), np.eye(m - 1)]), list(range(n - 1, n + m - 1))  # temp是一个对角矩阵，B是个递增序列
        mat = self.mat = np.hstack([self.mat, temp])  # 横向拼接
         #判断目标函数里是否还有负系数项
        while mat[0, 1:].min() < 0:  
            # 在目标函数里找到第一个负系数的变量，找到替入变量
            col = np.where(mat[0, 1:] < 0)[0][0] + 1  
            row = np.array([mat[i][0] / mat[i][col] if mat[i][col] > 0 else 0x7fffffff for i in
                            range(1, mat.shape[0])]).argmin() + 1  # 找到最严格约束的行，也就找到替出变量
            if mat[row][col] <= 0: return None  # 若无替出变量，原问题无界
            mat[row] /= mat[row][col]    #替入变量和替出变量互换
            ids = np.arange(mat.shape[0]) != row
            # 对i!= row的每一行约束条件，做替入和替出变量的替换
            mat[ids] -= mat[row] * mat[ids, col:col + 1]  
            B[row] = col  #用B数组记录替入的替入变量
        return mat[0][0] * (1 if self.max_mode else -1), {B[i]: mat[i, 0] for i in range(1, m) if B[i] < n} #返回目标值，对应x的解就是基本变量为对应的bi，非基本变量为0


"""
       minimize -x1 - 14x2 - 6x3
       st
        x1 + x2 + x3 <=4
        x1 <= 2
        x3 <= 3
        3x2 + x3 <= 6
        x1 ,x2 ,x3 >= 0
       answer :-32
    """
t = Simplex([-1, -14, -6])
t.add_constraint([1, 1, 1], 4)
t.add_constraint([1, 0, 0], 2)
t.add_constraint([0, 0, 1], 3)
t.add_constraint([0, 3, 1], 6)
print(t.solve())
print(t.mat)

调参侠求解LP

这里以为case2作为例子, 化为minimize标准型哦.

$min \ z = -x1 - x2 \\ s.t.$

$x1-x2<0 \\ -1.5x1 + x2 <= 0 \\ 50x1+20x2 <= 2000 \\ x1,x2 >=0$

api: from scipy.optimize import linprog

我一般都是用pycharm来看源码的api文档, 这样会印象深刻一下.

• c: 目标函数(minimize)的系数, 1D-array

• A_ub: 左边不等号相关的系数矩阵 2D-array

• b_ub: 即A_up中的对应行的右边值 1D-array

• A_eq: 左边等号相关的系数矩阵 2D-array

• b_eq: 即A_eq中的对应行的右边值 1D-array

• bounds: sequence

  • (min, max) pairs for each element in "X"

from scipy.optimize import linprog
import numpy as np


# 1. 目标函数系数
c = np.array([-1, -1])
# 2. 不等号
A_ub = np.array([[-1.5, 1], [50, 20]])
b_ub = np.array([0, 2000])
# 3. 等号(此题没有A_eq, b_eq)
# 4. bounds
limit_1, limit_2 = (0, None), (0, None)
# 5. solve: default "simplex"
ret = linprog(c, A_ub, b_ub, bounds(limit_1, limit_2))
print(ret)

     fun: -62.5
 message: 'Optimization terminated successfully.'
     nit: 2
   slack: array([0., 0.])
  status: 0
 success: True
       x: array([25. , 37.5])