📜  单纯形算法 - 表格法

📅  最后修改于: 2022-05-13 01:54:24.024000             🧑  作者: Mango

单纯形算法 - 表格法

单纯形算法是线性规划中著名的优化技术。
LPP(线性规划问题)的一般形式是
  Max/Min Z = c^tX  s.t. AX \leq b       X \geq 0

示例:让我们考虑以下最大化问题。
  Max x_1 + x_2  s.t.   x_1 + x_2 + x4 = 8  2x_1 + x_2 + x_3 = 10
初始构建步骤:

  • 建立你的矩阵 A。A 将包含约束的系数。
  • 矩阵 b 将包含资源量。
  • 矩阵 c 将包含目标函数或成本的系数。

针对上述问题——
矩阵 A – 在迭代 0

在迭代 0


表格说明——
B:基础,包含基本变量。单纯形算法从形成一个身份矩阵的那些变量开始。在上面例如 x4 和 x3 形成一个 2×2 单位矩阵。
CB:目标函数中基本变量的系数。目标函数不包含 x4 和 x3,因此它们为 0。
XB:资源数量或者我们可以说是RHS的约束。
yi :完整的矩阵 A。

单纯形算法1. 从与单位矩阵相关的初始基开始。 2. 计算相对利润。对于 MAX 问题 -如果所有相对利润都小于或等于 0,则当前基础是最优基础。停止。否则继续 3.对于 MIN 问题如果所有的相对利润都大于或等于 0,则当前基数是最优基数。停止。否则继续 3。 3. 找到对应于最大相对利润的列。假设列 k 具有最大 Rel。利润。所以 xk 将进入基础。 4. 执行最小比率测试以确定哪个变量将离开基础。    min ratio test:  XBr/y_{rk} = min\{XB_i/y_{ik}\}  最小元素的索引,即“r”将确定离开变量。索引 r 处的基本变量将离开基础。注意:最小比率测试总是对正元素进行。 5. 很明显,输入的变量不会形成单位矩阵,所以我们将不得不执行行操作使其再次成为单位。找到枢轴元素。索引 (r, k) 处的元素将是枢轴元素,而行 r 将是枢轴行。 6. 将第 r 行除以枢轴使其为 1。然后从其他行中减去 c*(rth row) 使它们为 0,其中 c 是使该行为 0 所需的系数。

迭代 1 中的表

迭代 1 中的表

相对利润的计算 – (Cj – Zj),其中 Cj 是 Z 中的系数,Zj 是 yi*CB
C1 – Z1 = 1 – (1*0 + 2*0)
C2 – Z2 = 1 – (1*0 + 1*0)
C3 – Z3 = 0 – (0*0 + 1*0)
C4 – Z4 = 0 – (1*0 + 0*0)

所以相对利润是- 1, 1, 0, 0(如表所示)
显然,并非所有的相对利润都小于或等于 0。因此将执行下一次迭代。
确定进入变量和离开变量。
指数 1 处的最大相对利润 1。因此 x1 将进入基础。
   min ratio test:  XBr/y_{rk} = min\{8/1, 10/2\}
(8, 5) 的最小值是 5,它在索引 2 处。所以 x3 将离开基础。

由于输入的 x1 执行所需的行操作来制作单位矩阵。

枢轴索引 = [2, 4]
枢轴元素 = 2

将第二行除以枢轴元素,即 2 使其为 1。
并从 R1 中减去 1*R2 使其为 0
见下表。

迭代 2 中的表

迭代 2 中的表

相对利润 = 0, 1/2, -1/2, 0
枢轴索引 = [1, 5]
枢轴元素 = 1/2
执行必要的行操作。
见下表

表在迭代 3


相对利润 = 0, 0, 0, -1
由于所有的相对利润都小于或等于 0。所以达到了最优。
这将是最终的单纯形表,也是最优的。
最优 Z 值 = 6*1 + 2*1 = 8

执行此算法时可能会发生以下情况。

  • 案例 1 – 无界解决方案
    如果对应于最大相对利润的列仅包含非正实数,那么我们将无法执行最小比率测试。因此,它被报告为无界解决方案。
  • 案例 2 – 替代解决方案
    如果在任何迭代中,非基本变量的相对利润中的任何一个结果为 0,则它包含替代解决方案。将存在许多最佳解决方案。

示例 2
上面的例子是一个相等的情况,我们能够找到初始的基础。现在我们将对一个没有恒等式形成的示例执行单纯形法。
  MAX 2x_1 + 5x_2  s.t. x_1 + x_2 \leq 6  x_2 \leq 3  x_1 + 2x_2 \leq 9
将上述问题转换为标准形式即
  MAX 2x_1 + 5x_2  s.t. x_1 + x_2 + x_3 = 6  x_2 + x_4 = 3  x_1 + 2x_2 + x_5 = 9
其中 x3、x4 和 x5 是松弛变量。这些将形成身份,从而形成最初的基础。
迭代 0 处的表

迭代 0 时的表

现在继续前面的例子。
迭代 1 中的表

迭代 1 中的表


相对利润 = 2, 5, 0, 0, 0
枢轴索引 = [2, 5]
枢轴元素 = 1

迭代 2 中的表

迭代 2 中的表


相对利润 = 2, 0, 0, -5, 0
枢轴索引 = [1, 4]
枢轴元素 = 1

迭代 3 中的表

迭代 3 中的表


相对利润 = 0, 0, 0, -2, -3, 0
因为所有的相对利润都小于等于 0。达到最优。
这是最终的单纯形表,也是最优的。
最优 Z 值 = 3*2 + 3*5 + 0*0 = 21

单纯形算法的代码实现

import numpy as np 
from fractions import Fraction # so that numbers are not displayed in decimal.
  
print("\n                 ****SiMplex Algorithm ****\n\n")
  
# inputs 
  
# A will contain the coefficients of the constraints
A = np.array([[1, 1, 0, 1], [2, 1, 1, 0]])
# b will contain the amount of resources 
b = np.array([8, 10])           
# c will contain coefficients of objective function Z      
c = np.array([1, 1, 0, 0])             
  
# B will contain the basic variables that make identity matrix
cb = np.array(c[3])
B = np.array([[3], [2]])          
 # cb contains their corresponding coefficients in Z   
cb = np.vstack((cb, c[2]))        
xb = np.transpose([b])                 
# combine matrices B and cb
table = np.hstack((B, cb))             
table = np.hstack((table, xb))         
# combine matrices B, cb and xb
# finally combine matrix A to form the complete simplex table
table = np.hstack((table, A))         
# change the type of table to float
table = np.array(table, dtype ='float') 
# inputs end
  
# if min problem, make this var 1
MIN = 0
  
print("Table at itr = 0")
print("B \tCB \tXB \ty1 \ty2 \ty3 \ty4")
for row in table:
    for el in row:
                # limit the denominator under 100
        print(Fraction(str(el)).limit_denominator(100), end ='\t') 
    print()
print()
print("Simplex Working....")
  
# when optimality reached it will be made 1
reached = 0     
itr = 1
unbounded = 0
alternate = 0
  
while reached == 0:
  
    print("Iteration: ", end =' ')
    print(itr)
    print("B \tCB \tXB \ty1 \ty2 \ty3 \ty4")
    for row in table:
        for el in row:
            print(Fraction(str(el)).limit_denominator(100), end ='\t')
        print()
  
    # calculate Relative profits-> cj - zj for non-basics
    i = 0
    rel_prof = []
    while i0:
            flag = 1
            break
        # if all relative profits <= 0
    if flag == 0:
        print("All profits are <= 0, optimality reached")
        reached = 1
        break
  
    # kth var will enter the basis
    k = rel_prof.index(max(rel_prof))
    min = 99999
    i = 0;
    r = -1
    # min ratio test (only positive values)
    while i0 and table[:, 3 + k][i]>0): 
            val = table[:, 2][i]/table[:, 3 + k][i]
            if val

对于以上内容,只需插入所需的值,您将通过单纯形算法获得 LPP 的详细逐步解决方案。