线性规划单纯形法python实现与代码详细解读
本文为大家详细解读线性规划单纯形法的原理,以及超详细的python源码解读!
线性规划单纯形法python实现与代码详细解读
1 单纯形法(Simplex method)
- 单纯形法是线性规划求解的经典方法之一,它的理论基础在于线性规划的解一定可以在顶点,即基可行解处取得。
- 举个简单例子,大家熟悉的二元线性规划中,目标函数的直线在坐标系上移动,在与可行域的交界处会取得最优解,即使目标函数与边界线相重合,也必然会有顶点位于该重合的直线上,即在顶点处必然会取得最优解。
- 基于以上论述,单纯形法就是从一个基可行解开始,通过一次改变一个基变量,从而实现从当前顶点转至下一个顶点的过程,结合基变量的模型表达式如下:
m a x c B X B + c N X N max\quad c_B X_B+c_N X_N maxcBXB+cNXN
s . t . B X B + N X N = b s.t.\quad B{X}_{B}+N{X}_{N}=b s.t.BXB+NXN=b
X B ≥ 0 , X N ≥ 0 \qquad X_{B}\ge0,X_{N}\ge0 XB≥0,XN≥0
2 编程思路
对于单纯形法的编程实现,其主要基于单纯形表,即约束系数矩阵,约束常数项以及检验数组成,通过初等行变换不断构建列线性无关的基,观察检验数的变化,当检验数表示其他变量入基不会带来更好的目标时,也就达到了最优解,具体的过程解读以及内在含义将在下方结合代码说明。
3 python实现原理解读
鉴于博主也为算法小白,源代码参考于运筹系列1:线性规划单纯形法python代码
原博主的代码实现方式非常简洁,但对于初学者而言还是会相对难以理解,本文也旨在对其做较为详细的说明与解读。
- 数据准备
该算法的输入为一个已有一组基的单纯形表,目标函数取最大值,如下表
( 1 14 6 0 0 0 0 0 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 26 ( x 1 x 2 x 3 x 4 x 5 x 6 x 7 b ) ) \quad\quad\quad\quad\quad\quad\quad\quad\begin{pmatrix} 1&14&6&0&0&0 &0&0\\\\ 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&26\\\\ (x_1&x_2&x_3&x_4&x_5&x_6 &x_7&b)\\\\ \end{pmatrix} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛11100(x1141003x261011x301000x400100x500010x600001x7042326b)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
在该矩阵中,第一行(row)的含义为目标函数值由非基变量表达的式子,即 Z = σ j X j + Z 0 Z=\sigma_jX_j+Z_0 Z=σjXj+Z0,其第一行即为 Z = X 1 + 14 X 2 + 6 X 3 + 0 Z=X_1+14X_2+6X_3+0 Z=X1+14X2+6X3+0 注意,第一行最后一列为 Z 0 Z_0 Z0,由于该算例的初始基可行解均为辅助变量,即 X 4 , X 5 , X 6 , X 7 X_4,X_5,X_6,X_7 X4,X5,X6,X7 ,则刚好原问题中的初始变量的 c j c_j cj都得以保留。(在一般问题当中也较为普遍,若每一行都添加了辅助变量来构造初始基可行解)
由于此时目标函数由一系列非基变量的线性组合构成,即与基变量的取值无关,则可通过该表达式观察是否有更好的解,在此例中,可观察 X 2 X_2 X2的系数最大,即 X 2 X_2 X2增加一个单位,目标函数值即可增加14个单位,则 X 2 X_2 X2应当入基从而获得更优基。 没错,这一行的实际意义就是检验数 σ j = c j − ∑ i = 1 m c m a i j \sigma_j =c_j-\sum_{i=1}^mc_ma_{ij} σj=cj−∑i=1mcmaij, 即该非基变量增加一个单位,其收益 c j c_j cj与在约束条件下,基变量随其增加1而相应减少的收益 ∑ i = 1 m c m a i j \sum_{i=1}^mc_ma_{ij} ∑i=1mcmaij之差,这个式子也好理解,比如看矩阵第五行,即第四个约束条件,若 X 2 X_2 X2增加一个单位,则为了满足约束,基变量 X 7 X_7 X7需减少3个单位,其收益损失 c 7 ∗ a 57 = 0 ∗ 3 = 0 c_7*a_{57}=0*3=0 c7∗a57=0∗3=0个单位,其余基变量也同理,这也就是检验数的实际效果。
根据以上推断,在编写代码时,我们要做到不断的变换基变量,直至:
- 第一行表达式中的非基变量如果入基,则其需要离开表达式,保证表达式中没有基变量,否则变换非基变量,基变量同样也会改变,则目标函数值变化方向不定,检验数失去效果。
- 第一行的前七个数均小于零,即满足最优性检验条件
有了这样的准则,我们只需编写数据读取与初始化操作,出基入基操作(pivot),与最优性判断条件(嵌在solve方法中)。
4 python代码
- 数据读取与初始化
数据如下:
1 14 6 0 0 0 0 0
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
##matrix即为上面展示的矩阵
matrix = np.loadtxt("data.txt", dtype=float)
(b_num, c_num) = matrix.shape
# b_num为约束行数加第一行的目标函数表达式行,c_num为约束系数列加上常数项b列
indexes_B = list(range(c_num - b_num, c_num - 1)) #基变量列表
solve()
printSol()
- 出基入基操作Pivot()
def pivot():
# 求转入基变量的列索引id
l = list(matrix[0][:-1])
index_intoB = l.index(max(l))#即取检验数的最大值列索引
# 算出基变量在表中的行索引,用b列除以入基变量列后的值进行比较判断
col_inB = []#入基的基变量的列,看哪个基变量出基
for i in range(b_num):
if matrix[i][index_intoB] == 0:
col_inB.append(0.)
else:
col_inB.append(matrix[i][-1] / matrix[i][index_intoB])
index_outB = col_inB.index(min([x for x in col_inB[1:] if x > 0]))#取b/aij的最小值,否则其他约束会不可行
#将基变量列表中进行新旧基的替换
indexes_B[index_outB - 1] = index_intoB #index_outB - 1为实际出基变量在基变量列表中的索引,因为矩阵第一行是c,约束矩阵A从第二行开始
cur_cell_to1 = matrix[index_outB][index_intoB] #新基的应为1的位置当前的数值
matrix[index_outB] /= cur_cell_to1 #该行全部除以该值,使其为1,即该行标准化
#通过行变换将其他行非零元去除
#注意,这一操作对第一行也执行了,是在将表达式中的新基变量也给消除!用其他非基变量来表示!第一行的最后一个数值极为z0,即目标函数值!
for row in [x for x in range(b_num) if x != index_outB]:
cur_cell_to0 = matrix[row][index_intoB]
matrix[row] -= cur_cell_to0 * matrix[index_outB]
- 算法执行与最优性检验
def solve():
flag = True
while flag:
if max(list(matrix[0][:-1])) <= 0: # 直至所有系数小于等于0
flag = False
else:
pivot()
def printSol():
for i in range(c_num - 1):
if i in indexes_B:
print("x" + str(i) + "=%.2f" % matrix[indexes_B.index(i) + 1][-1])
else:
print("x" + str(i) + "=0.00")
print("objective is %.2f" % (-matrix[0][-1]))
- 整体代码(考虑编译先后顺序)
import numpy as np
def pivot():
# 求转入基变量的列索引id
l = list(matrix[0][:-1])
index_intoB = l.index(max(l))#即取检验数的最大值列索引
# 算出基变量在表中的行索引,用b列除以入基变量列后的值进行比较判断
col_inB = []#入基的基变量的列,看哪个基变量出基
for i in range(b_num):
if matrix[i][index_intoB] == 0:
col_inB.append(0.)
else:
col_inB.append(matrix[i][-1] / matrix[i][index_intoB])
index_outB = col_inB.index(min([x for x in col_inB[1:] if x > 0]))#取b/aij的最小值,否则其他约束会不可行
#将基变量列表中进行新旧基的替换
indexes_B[index_outB - 1] = index_intoB #index_outB - 1为实际出基变量在基变量列表中的索引,因为矩阵第一行是c,约束矩阵A从第二行开始
cur_cell_to1 = matrix[index_outB][index_intoB] #新基的应为1的位置当前的数值
matrix[index_outB] /= cur_cell_to1 #该行全部除以该值,使其为1,即该行标准化
#通过行变换将其他行非零元去除
#注意,这一操作对第一行也执行了,是在将表达式中的新基变量也给消除!用其他非基变量来表示!第一行的最后一个数值极为z0,即目标函数值!
for row in [x for x in range(b_num) if x != index_outB]:
cur_cell_to0 = matrix[row][index_intoB]
matrix[row] -= cur_cell_to0 * matrix[index_outB]
def solve():
flag = True
while flag:
if max(list(matrix[0][:-1])) <= 0: # 直至所有系数小于等于0
flag = False
else:
pivot()
def printSol():
for i in range(c_num - 1):
if i in indexes_B:
print("x" + str(i) + "=%.2f" % matrix[indexes_B.index(i) + 1][-1])
else:
print("x" + str(i) + "=0.00")
print("objective is %.2f" % (-matrix[0][-1]))
matrix = np.loadtxt("data.txt", dtype=float)
(b_num, c_num) = matrix.shape
# b_num为约束行数加第一行的目标函数表达式行,c_num为约束系数列加上常数项b列
indexes_B = list(range(c_num - b_num, c_num - 1)) #基变量列表
solve()
printSol()
5 后记
博主研究方向是高速铁路运行图(大规模混合整数规划)的优化,现在正在啃整数规划的基础原理,这个单纯形法是我研学之路的开始,未来会对各类整数规划经典问题进行学习并分享我的经验,欢迎大家交流!
![Logo](https://devpress.csdnimg.cn/79de2bf0b7994defa4242ef90d5513fa.jpg)
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)