单纯形法和对偶单纯形法数学建模

检验数Ci - Zi

前言我在学习这个知识点的时候，花了不少时间，苦于网上的教程要么和老师描述的方法不一致（当然应该也可以），要么不全，所以当时在做作业题的时候，花了不少时间(>10hours)。因此想写一篇博文来清晰的记录这两种方法的用法。

原创链接: 单纯形法和对偶单纯形法, 转载请务必添加此声明!
https://blog.csdn.net/weixin_43850253/article/details/107793185

单纯形法 1. 停止条件
当检验数Ci - Zi 这一行的数都小于0时，可以停止。

2. 选择主元标准
我总结了一个小技巧: 先选检验数 Ci-Zi 那一行最大的。然后再选择 b / ai (ai > 0)最小的。
这里的ai的意思呢，大家看下面这个例子哈。

文章图片

我们先看 Ci - Zi 最大的那一列，是3，即 X2 那一列，这个时候我们知道换入主元是 X2。然后再看选择哪一行，这里X1 和 x5那两行可以考虑，因为呢(这里a1 = 2, a4 = 0, a5 = 4), a4 = 0不能考虑。b1 /a1 = 8 / 2 = 4, b5 / a5 = 12 / 4 = 3,我们选择小的那个。
因此换入主元是 x2, 换出主元是 x5 。

3. 一道例题
求m a xZ = 2 x 1 + 3 x 2 求~~~~max~Z = 2x_{1} + 3 x_{2} 求max Z=2x1?+3x2?
s . t . { x 1 + 2 x 2 ≤ 8 4 x 1 ≤ 16 4 x 2 ≤ 12 x i ≥ 0 , i = 1 , 2 s.t. \left\{ \begin{array}{l} x1 + 2x2 ≤ 8 \\ 4x1 ≤ 16 \\ 4x2 ≤ 12 \\ xi ≥ 0, i = 1,2 \end{array}\right. s.t.????????x1+2x2≤84x1≤164x2≤12xi≥0,i=1,2?
我们先对其标准化，即
求m a xZ = 2 x 1 + 3 x 2 求~~~~max~Z = 2x_{1} + 3 x_{2} 求max Z=2x1?+3x2?
s . t . { x 1 + 2 x 2 + x 3 = 8 4 x 1 + x 4 = 16 4 x 2 + x 5 = 12 x i ≥ 0 , i = 1 , 2 , 3 , 4 , 5 s.t. \left\{ \begin{array}{l} x1 + 2x2 + x3 = 8 \\ 4x1 + x4 = 16 \\ 4x2 + x5 = 12 \\ xi ≥ 0, i = 1,2,3,4,5 \end{array} \right. s.t.????????x1+2x2+x3=84x1+x4=164x2+x5=12xi≥0,i=1,2,3,4,5?

我们画出一个单纯形表，这里选择的主元的值为4(红色标记部分)。然后呢，这里基变量选择x3， x4， x5 是因为选他们方便，朋友们看一下基变量那一列就明白选他们计算量比较少。

检验数 Ci - Zi 有些大于0 。再计算。

检验数 Ci - Zi 有些大于0 。再计算。

	b	x1	x2	x3	x4	x5
x1	2	1	0	1	0	-1/2
x4	8	0	0	-4	1	2
x2	3	0	1	0	0	1/4
检验数Ci - Zi		0	0	-2	0	1/4

检验数 Ci - Zi 有些大于0 。再计算。

	b	x1	x2	x3	x4	x5
x1	4	1	0	0	1/4	0
x5	4	0	0	-2	1/2	1
x2	2	0	1	1/2	-1/8	0
检验数Ci - Zi		0	0	-3/2	-1/8	0

检验数 Ci - Zi 都小于0，停止计算。得出最优解 x1 = 4, x2 = 2, x5 = 4,
max Z = 2x1 + 3 x2 = 2 × 4 + 3 × 2 = 8 + 6 = 14

对偶单纯形法 1. 使用要求/条件
要求 b 那一列至少有一个数小于 0，检验数Ci - Zi 都小于0。

2. 选择主元标准
先看b那一列，挑个最小的，然后再看 Ci - Zi 那一行，选 Zi / ai (ai < 0)最小。
大家看下面这个例子哈。

文章图片
这里选择-2/7的理由: 先看b那一列最小的是 -6 / 7，即换出主元是 x6 确定了，然后我们看检验数那一列， Z3 / a3 = (-5/7) / (-1/7) = 5, Z5 / a5 = (-3/7) / (-2/7) = 3 / 2, 这里选小的，故换入主元为 x 5。

3. 使用流程
先用单纯形法运算到满足对偶单纯形法的使用条件后再用对偶单纯形法。

4. 一道例题和源代码

文章图片

解题过程:
(1)先加入松弛变量
{ 3 x 1 ? 2 x 2 + x 3 = 3 5 x 1 + 4 x 2 ? x 4 = 10 2 x 1 + x 2 + x 5 = 5 \left\{ \begin{array}{l} 3x1 - 2x2 + x3 = 3 \\ 5x1 + 4x2 - x4 = 10 \\ 2x1 + x2 + x5 = 5 \end{array} \right. ????3x1?2x2+x3=35x1+4x2?x4=102x1+x2+x5=5?

(2)初始的单纯形表:选定x3, x4,x5 作为基变量,这里这样选择会使得这一步计算比较少
（表1）

(3)这里用单纯形法换入换出变量,这里选择换出变量为x4，换入变量为x2, 因为b/ai最小(ai>0)
这一步选择22/3
（表2）

	b	x1	x2	x3	x4	x5
x1	1	1	-2/3	1/3	0	0
x4	5	0	22/3	-5/3	-1	0
x5	3	0	7/3	-2/3	0	1
检验数	空	0	1	-1	0	0

(4)上一步检验数不全小于0，继续单纯形。这一步选择7/22.
（表3）

	b	x1	x2	x3	x4	x5
x1	16/11	1	0	2/11	-1/11	0
x2	15/22	0	1	-5/22	-3/22	0
x5	31/22	0	0	-3/22	7/22	1
检验数	空	0	0	-17/22	3/22	0

换入换出后（表4）

	b	x1	x2	x3	x4	x5
x1	13/7	1	0	1/7	0	2/7
x2	9/7	0	1	-2/7	0	3/7
x4	31/7	0	0	-3/7	1	22/7
检验数	空	0	0	-5/7	0	-3/7

可以看到此时解仍然部位整数，不过此时检验数已经都小于0.下面就可以使用对偶单纯形法了。
对于x1 = 13/7, 找割平面 1/7 x3 + 2/7 x5 ≥ 6/7,
即 -1 / 7 x3 - 2 / 7x5 + x6 = - 6 / 7.

添加割平面，得
（表5）

	b	x1	x2	x3	x4	x5	x6
x1	13/7	1	0	1/7	0	2/7	0
x2	9/7	0	1	-2/7	0	3/7	0
x4	31/7	0	0	-3/7	1	22/7	0
x6	-6/7	0	0	-1/7	0	-2/7	1
检验数	空	0	0	-5/7	0	-3/7	0

然后选定主元-2/7.

(5)
（表6）

	b	x1	x2	x3	x4	x5	x6
x1	1	1	0	0	0	0	1
x2	0	0	1	-1/2	0	0	3/2
x4	-5	0	0	-2	1	0	11
x5	3	0	0	1/2	0	1	-7/2
检验数	空	0	0	-1/2	0	0	-3/2

（6）得出
（表7）

	b	x1	x2	x3	x4	x5	x6
x1	1	1	0	0	0	0	1
x2	5/4	0	1	0	-1/4	0	-5/4
x3	5/2	0	0	1	-1/2	0	-11/2
x5	7/4	0	0	0	1/4	1	-3/4
检验数	空	0	0	0	-1/4	0	-17/4

最优解为(1, 5/4, 5/2, 0, 7/4, 0)
很不幸，这里仍然不为整数,则对x5 = 7/4, 有1/4 x4 + 1/4 x6 ≥ 3/4,
即将 -1/4 x4 - 1/4 x6 + x7 = -3/4加入割平面.

（表8）

	b	x1	x2	x3	x4	x5	x6	x7
x1	1	1	0	0	0	0	1	0
x2	5/4	0	1	0	-1/4	0	-5/4	0
x3	5/2	0	0	1	-1/2	0	-11/2	0
x5	7/4	0	0	0	1/4	1	-3/4	0
x7	-3/4	0	0	0	-1/4	0	-1/4	1
检验数	空	0	0	0	-1/4	0	-17/4	0

(7) 化简得
（表9）

	b	x1	x2	x3	x4	x5	x6	x7
x1	1	1	0	0	0	0	1	0
x2	2	0	1	0	0	0	-1	-2
x3	4	0	0	1	0	0	-5	-2
x5	1	0	0	0	0	1	-1	1
x7	3	0	0	0	1	0	1	-4
检验数	空	0	0	0	0	0	-4	-1

最优解为(1, 2, 4, 1, 3, 0, 0)
因此maxZ 的最优解为 3x1 - x2 = 3 * 1 - 2 = 1

代码来自于我的一个同学，特别感谢他:（操作都封装成函数了，大家解题部分只需要看主函数即main那部分），这里用到的numpy操作，如果没学过numpy的可以瞅瞅官网文档, 或者numpy中文网文档

""" 对偶单纯形法，主元选择过程靠人工，由于对偶单纯形和单纯形只是选择主元不同，所以都可以适用 """ from collections.abc import Iterable from fractions import Fraction import numpy as np# 让numpy输出的Fraction是居中的str而不是repr Fraction.__repr__ = lambda self: str(self).center(6)table: np.ndarray var: list count = 0# 创建分数类型的numpy数组 def farray(li): return np.array([Fraction(*each) if isinstance(each, Iterable) else Fraction(each) for each in li])# 交换基变量和非基变量 def swap(row, col): global table global var var[row] = col table[row] /= table[row][col] for i in range(len(table)): if i == row: continue table[i] -= table[i][col] * table[row]# 加入新约束 def new_constraint(li): global table global var res = np.expand_dims(farray([0] * table.shape[0]), axis=1) table = np.hstack((table, res)) table = np.vstack((table, farray(li))) table[[-1, -2], :] = table[[-2, -1], :] var.append(table.shape[1] - 1)# 显示矩阵 def show(): global count head = ['', 'b'] head.extend([f'x{i}' for i in range(1, table.shape[1])]) max_len = max(map(lambda t: len(t), head)) li2 = [] for row in table: li2.append([str(col) for col in row]) max_len = max(max_len, max(map(lambda t: len(t), li2[-1]))) for i in range(len(li2) - 1): li2[i].insert(0, f'x{var[i]}') li2[-1][0] = 'check_number' max_len = max(7, max_len, len(li2[-1][0]) // 2) head = list(map(lambda t: t.center(max_len), head)) for i in range(len(li2)): li2[i] = list(map(lambda t: t.center(max_len), li2[i])) li2[-1][0] = li2[-1][0].center(max_len * 2 + 1) count += 1 print(f'{count})') print(*head) for row in li2: print(*row) print()if __name__ == '__main__': table = np.array([ farray([3, 3, -2, 1, 0, 0]), farray([10, 5, 4, 0, -1, 0]), farray([5, 2, 1, 0, 0, 1]), farray([0, 3, -1, 0, 0, 0]) ]) # 基变量 var = [3, 4, 5]show()# 第0行和第1列交换 swap(0, 1) show()swap(1, 2) show()swap(2, 4) show()new_constraint([(-6, 7), 0, 0, (-1, 7), 0, (-2, 7), 1]) show()swap(3, 5) show()swap(2, 3) show()new_constraint([(-3, 4), 0, 0, 0, (-1, 4), 0, (-1, 4), 1]) show()swap(4, 4) show()

【单纯形法和对偶单纯形法】下面是输出结果:

1) bx1x2x3x4x5 x31370102 x41030014 x2521001 check_number50001 2) bx1x2x3x4x5 x113/7101/702/7 x431/700-3/7122/7 x29/701-2/703/7 check_number00-5/70-3/7 3) bx1x2x3x4x5x6 x113/7101/702/70 x431/700-3/7122/70 x29/701-2/703/70 x6-3/700-4/70-1/71 check_number00-5/70-3/70 4) bx1x2x3x4x5x6 x17/410001/41/4 x419/4000113/4-3/4 x23/201001/2-1/2 x33/400101/4-7/4 check_number0000-1/4-5/4 5) bx1x2x3x4x5x6x7 x17/410001/41/40 x419/4000113/4-3/40 x23/201001/2-1/20 x33/400101/4-7/40 x7-3/40000-1/4-1/41 check_number0000-1/4-5/40 6) bx1x2x3x4x5x6x7 x111000001 x4-500010-413 x2001000-12 x3000100-21 x53000011-4 check_number00000-1-1 7) bx1x2x3x4x5x6x7 x111000001 x65/4000-1/401-13/4 x25/4010-1/400-5/4 x35/2001-1/200-11/2 x57/40001/410-3/4 check_number000-1/400-17/4 8) bx1x2x3x4x5x6x7x8 x1110000010 x65/4000-1/401-13/40 x25/4010-1/400-5/40 x35/2001-1/200-11/20 x57/40001/410-3/40 x8-1/2000-1/200-1/21 check_number000-1/400-17/40 9) bx1x2x3x4x5x6x7x8 x1110000010 x63/2000001-3-1/2 x23/2010000-1-1/2 x33001000-5-1 x53/2000010-11/2 x410001001-2 check_number000000-4-1/2 10) bx1x2x3x4x5x6x7x8x9 x11100000100 x63/2000001-3-1/20 x23/2010000-1-1/20 x33001000-5-10 x53/2000010-11/20 x410001001-20 x9-1/20000000-1/21 check_number000000-4-1/20 11) bx1x2x3x4x5x6x7x8x9 x11100000100 x62000001-30-1 x22010000-10-1 x34001000-50-2 x51000010-101 x4300010010-4 x8100000001-2 check_number000000-40-1