线性规划的单纯形表—R和Pyhton实现
线性规划是运筹学中的一个基本分支，它广泛应用现有的科学技术和数学方法，解决实际中的问题，帮助决策人员选择最优方针和决策，自1947年丹捷格提出了一般线性规划问题求解的方法———单纯形法之后，线性规划在理论上趋向成熟，特别是在电子计算机能处理成千上万个约束条件和决策的线性规划问题之后，线性规划的适用领域更为广泛，是对有限的资源进行合理分配，企业提高生产效率，从而获得最佳经济效益的有效工具。
一、问题的提出
例：（生产计划问题）假设某厂计划生产甲、乙两种产品，其主要原材料有钢材360kg，铜材300kg及专用设备能力200台时，已原材料和设备的单间消耗定额以及单位产品所获利润如下表所示：
| 甲 | 乙 | 现有资源 | |
|---|---|---|---|
| 钢材 | 9 | 4 | 360 |
| 铜材 | 3 | 10 | 300 |
| 设备台时 | 4 | 5 | 200 |
| 利润 | 60 | 120 |
问如何安排生产方使该厂所获利润最大？
二、线性规划模型
设生产甲乙两种产品的数量为\(x_1\)和\(x_2\)，建立数学模型如下：
\[max \ \ \ z=60x_1+120x_2\\\ \ \ \begin{cases} 9x_1+4x_2\leq360\\3x_1+10x_2\leq300\\4x_1+5x_2\leq200\\ x_1\geq0;x_2\geq0 \end{cases}
\]
标准化：
\[max \ \ \ z=60x_1+120x_2\\\ \ \ \begin{cases} 9x_1+4x_2+x_3=360\\3x_1+10x_2+x_4=300\\4x_1+5x_2+x_5=200\\ x_1\geq0;x_2\geq0 \end{cases}
\]
三、R计算程序
3.1 R程序
Lp <-function(c,A,b,lav)
#c为目标函数系数向量；A为添加松弛变量后的系数矩阵；
#b为常向量;lav松弛变量序号（构成单位阵）
{
##求初始基可行解、检验数
n = ncol(A);m = nrow(A)
sol = rep(0,n)
for(i in 1:length(lav))sol[lav[i]] = b[i]
B = lav #基变量向量序号（若1,3为基变量序号，则B=c(1,3))
sigma = c - c[lav]%*%A
print("初始单纯形表！")
DF = data.frame("CB" = c[B],"Base" = B,"b"=b, "x" = A)
print(list("frame" = DF,"sigma" = sigma))
while(any(sigma>0))
{
infty = which(sigma>0)
for(i in infty)if(all(A[,i] <= 0))print("答：存在无界解！")
psigma = max(sigma)
pcol = which(sigma == psigma) #主元素所在列
sita0 = b/A[,pcol]
sita = min(sita0)
prow = which(sita0 == sita) #主元素所在行
B[prow] = pcol #换基变量
##列出新单纯形表
b[prow] = b[prow]/A[prow,pcol] #对主元素行
A[prow,] = A[prow,]/A[prow,pcol]
for(i in c(1:m)[-prow])
{
b[i] = b[i] - b[prow]*A[i,pcol]
A[i,] = A[i,] - A[prow,]*A[i,pcol]
}
sigma = c - c[B]%*%A
print("过程单纯形表！")
DF = data.frame("CB" = c[B],"Base" = B,"b"=b, "x" = A)
print(list("frame" = DF,"sigma" = sigma))
}
if(any(sigma[-B]==0))
{print("答：存在无穷多最优解！")
}else{print("答：唯一最优解！")}
print("最终单纯形表！")
DF = data.frame("CB" = c[B],"Base" = B,"b"=b, "x" = A)
print(list("goal coef" = c,"frame" = DF,"sigma" = sigma))
sol = rep(0,n)
for(i in 1:length(B))sol[B[i]] = b[i]
return(list("最优解"=sol,"最大值"=sum(sol*c)))
}
#输入数据
c=c(60,120,0,0,0)
A=matrix(c(9,3,4,4,10,5,1,0,0,0,1,0,0,0,1),3)
b=c(360,300,200)
lav=c(3,4,5)
#运行结果
Lp(c,A,b,lav)
3.2 计算步骤展现
#初始单纯形表
$frame
CB Base b x.1 x.2 x.3 x.4 x.5
1 0 3 360 9 4 1 0 0
2 0 4 300 3 10 0 1 0
3 0 5 200 4 5 0 0 1
$sigma
[,1] [,2] [,3] [,4] [,5]
[1,] 60 120 0 0 0`
#过程单纯形表1
$frame
CB Base b x.1 x.2 x.3 x.4 x.5
1 0 3 240 7.8 0 1 -0.4 0
2 120 2 30 0.3 1 0 0.1 0
3 0 5 50 2.5 0 0 -0.5 1
$sigma
[,1] [,2] [,3] [,4] [,5]
[1,] 24 0 0 -12 0
#过程单纯形表2
$frame
CB Base b x.1 x.2 x.3 x.4 x.5
1 0 3 84 0 0 1 1.16 -3.12
2 120 2 24 0 1 0 0.16 -0.12
3 60 1 20 1 0 0 -0.20 0.40
$sigma
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 -7.2 -9.6
#最终单纯形表
$frame
CB Base b x.1 x.2 x.3 x.4 x.5
1 0 3 84 0 0 1 1.16 -3.12
2 120 2 24 0 1 0 0.16 -0.12
3 60 1 20 1 0 0 -0.20 0.40
$sigma
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 -7.2 -9.6
$最优解
[1] 20 24 84 0 0
$最优值
[1] 4080
四、Python计算程序
4.1 Pyhton程序
import pandas as pd
import numpy as np
"""
系数矩阵的形式:
b x1 x2 x3 x4 x5
obj 0 70 30 0 0 0
x3 540 3 9 1 0 0
x4 450 5 5 0 1 0
x5 720 9 3 0 0 1
①第一行是目标函数的系数；亦即各变量对应的检验数；第2~4行是约束条件的系数
②第一列是约束方程的常数项
③对于目标函数的更新我们同样采用矩阵的变换,所以obj对应的第一列表示的是目标函数的相反数
"""
"""
运行如下代码后得到的结果如下所示
最终的最优单纯性法是:
b x1 x2 x3 x4 x5
obj -5700 0 0.0 0 -2.0 -6.666667
x3 180 0 0.0 1 -2.4 1.000000
x2 15 0 1.0 0 0.3 -0.166667
x1 75 1 0.0 0 -0.1 0.166667
目标函数的值: 5700
最优决策变量是:
x1 = 75
x2 = 15
"""
matrix = pd.DataFrame(
data=np.array([
[0, 60, 120, 0, 0, 0],
[360, 9, 4, 1, 0, 0],
[300, 3, 10, 0, 1, 0],
[200, 4, 5, 0, 0, 1]
]),
index=['obj', 'x3', 'x4', 'x5'],
columns=['b', 'x1', 'x2', 'x3', 'x4', 'x5']
)
i=1
print("第",i,"个单纯性表是:")
print(matrix)
# 判断检验数是否大于0
c = matrix.iloc[0, 1:]
while c.max() > 0:
c = matrix.iloc[0, 1:]
i=i+1
# 选择入基的变量
in_x = c.idxmax() # 该函数运行的结果是x1,即入基变量
in_x_v = c[in_x] # 得到入基变量的系数
# 选择出基变量,即要计算θ的值,并比较大小
b = matrix.iloc[1:, 0]
# 选择入基变量对应的列
in_x_a = matrix.iloc[1:][in_x]
# 得到出基变量
out_x = (b / in_x_a).idxmin()
# 完成入基和出基的操作,即对矩阵做初等行变化
matrix.loc[out_x, :] = matrix.loc[out_x, :] / matrix.loc[out_x][in_x]
for idx in matrix.index:
if idx != out_x:
matrix.loc[idx, :] = matrix.loc[idx, :] - matrix.loc[out_x, :] * matrix.loc[idx, in_x]
# 索引变换
index = matrix.index.tolist() # 得到所以行索引的值
i = index.index(out_x) # 得到出基变量的下标
index[i] = in_x
matrix.index = index
print("第",i,"个单纯性表是:")
print(matrix)
# 输出结果
print("最终的最优单纯性法是:")
print(matrix)
print("目标函数的值:", - matrix.iloc[0, 0])
# 列的数量减去行的数量得到非基变量的数量,且非基变量一定是下标在前的变量
# 比如非基变量的个数为2,则非基变量一定是x1,x2
print("最优决策变量是:")
x_count = (matrix.shape[1] - 1) - (matrix.shape[0] - 1)
X = matrix.iloc[0, 1:].index.tolist()[:x_count]
for xi in X:
print(xi, "=", matrix.loc[xi, 'b'])
4.2 计算步骤展现
第 1 个单纯性表是:
b x1 x2 x3 x4 x5
obj 0 60 120 0 0 0
x3 360 9 4 1 0 0
x4 300 3 10 0 1 0
x5 200 4 5 0 0 1
第 2 个单纯性表是:
b x1 x2 x3 x4 x5
obj -3600 24.0 0 0 -12.0 0
x3 240 7.8 0 1 -0.4 0
x2 30 0.3 1 0 0.1 0
x5 50 2.5 0 0 -0.5 1
第 3 个单纯性表是:
b x1 x2 x3 x4 x5
obj -4080 0.0 0 0 -7.20 -9.60
x3 84 0.0 0 1 1.16 -3.12
x2 24 0.0 1 0 0.16 -0.12
x1 20 1.0 0 0 -0.20 0.40
第 3 个单纯性表是:
b x1 x2 x3 x4 x5
obj -4080 0.0 0 0 -7.20 -9.60
x3 84 0.0 0 1 1.16 -3.12
x2 24 0.0 1 0 0.16 -0.12
x1 20 1.0 0 0 -0.20 0.40
最终的最优单纯性法是:
b x1 x2 x3 x4 x5
obj -4080 0.0 0 0 -7.20 -9.60
x3 84 0.0 0 1 1.16 -3.12
x2 24 0.0 1 0 0.16 -0.12
x1 20 1.0 0 0 -0.20 0.40
目标函数的值: 4080
最优决策变量是:
x1 = 20
x2 = 24
五、总结
运用所学运筹学知识，针对该公司生产计划提出一些科学决策方案，从而达到资源充分利用的目的。通过对方案的提出、分析和解决对策的制定，使我们能够运用运筹学知识和相关工具解决一些实际性问题，加深对该课程的认识。同时，通过软件的计算使用，能够实现将理论与实践相结合的目的，增强我们动手操作能力和工作协调力。
