首先,讲讲思路
我们之前通过求基解的方式得到了可行域的极点(可用的极点),而且我们知道在我们所有求到的极点中必定有一个最优解,我们只需要逐个比较就可以了。
但有些时候,时间开销会非常大,这种方法并不一定合适。
我们于是可以换个思路:
若从某一基本可行解出发,每次总是寻求比上一个更“好”的基本可行解,直至找到最优解。这样就可以大大减少计算量,其实这样的思想应用在非常多的地方,比如梯度下降等等。
于是,我们需要解决下面这三个问题:
(1) 如何判断当前的基本可行解是最优解(迭代终点);
(2) 如何寻找下一个改善的基本可行解(迭代关系式);
(3) 如何得到一个初始的基本可行解(迭代起点);
于是,就有了我们的单纯形方法:
基本过程:
第一步:构造一个初始的基本可行解。
第二步:判断当前基本可行解是否为最优解。
第三步:若当前解不是最优解,则要进行基变换迭代到下一个基本可行解。
所以,我们可以写出如下代码:
import numpy as np
if __name__=="__main__":
n=int(input('请输入有多少个变量'))
m=int(input('请输入有多少个表达式'))
print('请输入约束矩阵A')
A=[]
B=[]
C=[]
for i in range(m):
A.append([])
a=input().split()
for j in range(n):
A[i].append(int(a[j]))
a=input('请输入常系数\n').split()
for i in range(m):
B.append(int(a[i]))
a=input('请输入待求方程式\n').split()
C.append(0)
for i in range(n):
C.append(int(a[i]))
B0=[]
D=[]
for i in range(m):
D.append(0)
B0.append([])
for j in range(m):
if i==j:
B0[i].append(1)
else:
B0[i].append(0)
a=np.array(np.mat(A))
b=np.array(np.mat(B).T)
b0=np.array(np.mat(B0))
c=np.array(np.mat(C))
d=np.array(np.mat(D))
s=np.hstack([b,a])
s=np.vstack([c,s])
ss=np.vstack([d,b0])
s=np.hstack([s,ss])
s=np.array(s,dtype='float')
while s[0,1:].min()<0:
col=np.where(s[0,1:]<0)[0][0]+1
row=np.array([s[i][0]/s[i][col]if s[i][col]>0 else 0x7fffffff for i in
range(1,s.shape[0])]).argmin()+1
if s[row][col]<=0:print('无解')
s[row] /= s[row][col]
ids=np.arange(s.shape[0])!=row
s[ids]-=s[row]*s[ids,col:col+1]
# print(row)
print(s)
print('最小值为'+str(float(s[0][0]*-1)))
能够找到最小值,洒洒水啦,各位看看就好,我还需要进行优化修改