运筹学实验_单纯形法

实验目的：

python实现单纯形法、lingo实现灵敏度分析、对偶变量、人工变量

实验要求：

报告内容：

python实现：

代码：

import numpy as np
from fractions import Fraction as f 
#本程序仅适用于标准形式的线性规划问题
#其中目标函数求最大，约束条件等式化，决策变量非负
#本程序实现的是单纯形法和大M法
def judge(matrix):#最优性检验
    if max(matrix[-1][:-1]) <= 0:  
        flag = False
    else:
        flag = True
    return flag

def initialTable(a,matrix,vect,c):  #基变换
    c2=np.array(c)
    for i in range(len(a[-1])-1):
        matrix[-1][i]=c2[i]-np.dot(matrix[:-1,i],[c2[vect[j]-1]for j in range(len(vect))])
    a = [list(i) for i in matrix]
    return a,matrix

def simplexTable(matrix,vect,m,n):   #输出单纯形表
    print('*'*20)
    print(' ',end='\t')
    for i in range(n):
        print('x{}'.format(i+1),end='\t')
    print('b')
    for i in range(m+1):
        if i <= m-1:
            print('x{}'.format(vect[i]),end='\t')
        elif i == m:
            print('cj-zj',end='\t')
        for j in matrix[i]:
            print(f(str(j)).limit_denominator(),end='\t')   #输出分数形式
        print(end='\n')
        
def trans(a,matrix,vect,c):  #基变换
    #计算检验数
    maxSigmaj = max(matrix[-1][:-1]) #最大检验数
    k = a[-1].index(maxSigmaj)   #入基变量的下标
    l = {}                       #出基变量字典
    for i in a[:-1]:   
        if i[k] >0: #第k列第i行的值
            l[i[-1]/i[k]] = a.index(i)
    if len(l)!=0:
        theta = l[min(l)]    #出基变量的下标
    
    matrix[theta] = matrix[theta]/matrix[theta][k]
    for i in range(len(a)):
        if i != theta:
            matrix[i] = matrix[i] - matrix[i][k]*matrix[theta]    
vect[theta] = k+1   #基变量下标变化
    #计算z
    c2=[0 for i in range(len(c))]
    c3=np.array(c)
    for i in range(len(c)):#滤去-M
        if c[i]!=-M:
            c2[i]=c[i]
    c2=np.array(c2)
    x=[0 for i in range(len(c))]
    for i in range(len(vect)):
        x[vect[i]-1]=matrix[i][-1]
    x=np.array(x)
    matrix[-1][-1]=np.dot(c2,x)
    #把原来的列表也同时变换掉，为了方便索引
    a = [list(i) for i in matrix] 
    return a,matrix,vect

def print_solution(matrix,vect,c,m,n):#打印最优解
    print('*'*20)
    print('\n')
    print("最优解为：")
    for i in range(m):
        print('x{}*={}'.format(vect[i],matrix[i][-1]),end='，')
    print('\nz*={}'.format(matrix[-1][-1]))
    
def common(a,b,c,m,n):
    a=[i+[j] for i,j in zip(a,b)]+[c+[0]]
    matrix = np.array(a, dtype=np.float64)   
    vect = [int(input('输入基变量下标')) for i in range(m)]
    a,matrix=initialTable(a,matrix,vect,c)
    #输出单纯形表
    simplexTable(matrix,vect,m,n)
    while judge(matrix):
        a,matrix,vect = trans(a,matrix,vect,c)
        simplexTable(matrix,vect,m,n)
    print_solution(matrix,vect,c,m,n) 
def test():
    m = 4            #约束条件数
    n = 6             #决策变量数
    a = [[2,2,1,0,0,0],[1,2,0,1,0,0],[4,0,0,0,1,0],[0,4,0,0,0,1]]#资源矩阵
    b = [12,8,16,12]    #右端项
    c = [2,3,0,0,0,0]     #价值向量
    common(a,b,c,m,n)
def test1():    
    print("第一章第五节ppt例6 \n")
    m1 = 3#约束条件数
    n1 = 7          #决策变量数
    a1 = [[1,1,1,1,0,0,0],[-2,1,-1,0,-1,1,0],[0,3,1,0,0,0,1]]#资源矩阵
    b1 = [4,1,9] #右端项
    c1 = [-3,0,1,0,0,-M,-M] #价值向量
    common(a1,b1,c1,m1,n1)
def test2():  
    print("第一章习题7 \n")
    #这一题是无界解，但我们不做这方面的检查，程序会报错
    m2 = 3     #约束条件数
    n2 =  8     #决策变量数
    a2 = [[1,1,1,-1,1,0,0,0],[-2,0,1,0,0,-1,1,0],[0,-2,-1,0,0,0,0,1]]#资源矩阵
    b2 = [6,2,0]    #右端项
    c2= [2,-1,2,0,-M,0,-M,0]                 #价值向量
    common(a2,b2,c2,m2,n2)
def test3():
    print("第一章习题7 \n")
    #这一题是无界解，但我们不做这方面的检查，程序会报错
    m2 = 3                           #约束条件数
    n2 =  8                           #决策变量数
    a2 = [[1,1,1,-1,1,0,0,0],[-2,0,1,0,0,-1,1,0],[0,-2,-1,0,0,0,0,1]]#资源矩阵
    b2 = [6,2,0]                        #右端项
    c2= [2,-1,2,0,-M,0,-M,0]                    #价值向量
    common(a2,b2,c2,m2,n2)
def test4():
    print("第一章ppt 例子11 \n")
    m=3
    n=9
    a=[[2,1,2,0,0,3,1,0,0],[0,2,1,0,1,0,0,1,0],[1,0,0,2,1,0,0,0,1]]
    b=[200,300,200]
    c=[-1,-1,-1,-1,-1,-1,-M,-M,-M]
    common(a,b,c,m,n)
def main():
    #global m,n    
    global M
    M=1000
    test()
main()

Lingo实现：

代码：

model:
sets:
row/1..4/:b;
col/1..2/:x;
value/1..2/:c;
links(row,col):a;
endsets
data:
a=2 2
  1 2
  4 0
  0 4;
b=12 8 16 12;
c=2 3;
enddata
max=@sum(value(i):c(i) * x(i));
@for(row(i):@sum(col(j):a(i,j) * x(j))<=b(i));
@for(col(j):x(j)>=0);
end