对偶单纯形方法python实现

1.简介

通过采用表格形式的对偶单纯形方法类似的解题步骤，用python编程实现，该程序以最小化问题的求解过程为方法，做到输入约束条件以及初始判别数（当目标函数为最小化问题时，输入值为目标函数系数的相反数；当目标函数为最大化问题时，输入值为目标函数系数），输出单纯形表和最优解以及最优值，以及对偶问题的最优解。

2.程序设计思路

表格形单纯形法（以极小化为例）的步骤：

1. 给定一个初始对偶可行的基本解
2. 若b=B-1b≥0则停止计算，现行对偶可行的基本解就是最优解，否则，令br=minibi<0
3. 若对所有j，yrj≥0,则停止计算，原问题无可行解，否则，令zk-ckyrk=minjzj-cjyrj|yrj<0
4. 以yrj为主元，进行主元消去，返回步骤2）

 

3.程序代码

 

'''
对偶单纯形法
'''
from fractions import Fraction as f

import numpy as np

def getinput():
    global m,n,h
    string = input('''
    输入初始单纯形表形如
    3 2 1 0 18;-1 4 0 1 8;-2 1 0 0 0
    前m行表示m个约束的增广矩阵，最后一行表示判别数（
    当目标函数为最小化问题时，输入值为目标函数系数的相反数；
    当目标函数为最大化问题时，输入值为目标函数系数）
    输入：''')
    a = [list(map(eval,row.split())) for row in string.split(';')]
    matrix = np.mat(a)
    m,n = matrix.shape
    h=input('''
    若目标函数为极小化问题，则输入0；极大化问题，输入1
    ''')
    x=[]
    if h =='1':
        print('\n\n输入的目标函数为')
        for j in range(n):
            if matrix[-1, j]!=0:
                x.append(f'{matrix[-1, j]}*x_{j+1}')
        print('max z = '+' + '.join(x))
    else:
        print('\n\n输入的目标函数为')
        for j in range(n):
            if matrix[-1, j] != 0:
                x.append(f'{-matrix[-1, j]}*x_{j + 1}')
        print('min z = ' + ' + '.join(x))
    print('\n\n输入的方程为')
    for i in range(m-1):
        x = [f'{matrix[i,j]}*x_{j+1}' for j in range(n)]
        print(' + '.join(x),f'={matrix[i,-1]}')
    return matrix

def judge(matrix):
    mini = matrix[0][-1]
    index=0
    for i in range(m):
        if matrix[i][-1] < mini:
            mini = matrix[i][-1]
            index = i
    if matrix[index][-1] >= 0:  # 最后一列除了最后一行
        flag = False
    else:
        flag = True
    return flag
def pr(matrix,vect):   #输出单纯形表
    print('*'*20)
    print(' ',end='\t')
    for i in range(n-1):
        print('X_{}'.format(i+1),end='\t')
    print('b')
    for i in range(m):
        if i <= m-2:
            print('x_{}'.format(vect[i]),end='\t')
        elif i == m-1:
            print('z-c',end='\t')
        #for j in matrix[i]:
        for j in range(n):
            print(f(str(matrix[i][j])).limit_denominator(),end='\t')   #输出分数形式
        print(end='\n')
def trans(a,matrix,vect):  #转轴
    mini =matrix[0][-1]
    index=0
    # print(mini)
    for i in range(m-1):
        if matrix[i][-1]<mini:
            mini=matrix[i][-1]
            index=i
    l = {}
    #index=index-1
    for j in range(n-1):
        if matrix[index][j]<0:
            l[matrix[-1][j]/matrix[index][j]]=j
    pivot=l[min(l)]#入基变量的足标j(列）
    vect[index] = pivot + 1
    matrix[index] = matrix[index] / matrix[index][pivot]
    for i in range(m):
        if i != index:
            matrix[i]=matrix[i]-matrix[i][pivot]*matrix[index]
    a = [list(i) for i in matrix]  #把原来的列表也同时变换掉，为了方便索引
    return a,matrix,vect
def print_solution(matrix,vect):
    print('*'*20)
    print("原问题的最优解：")
    for i in range(1,n):
        if i in vect:
            print('x{}*={}'.format(i,f(str(matrix[vect.index(i)][-1])).limit_denominator()),end='，')
        else:
            print('x{}*={}'.format(i,0),end='，')
    print('\n')
    print("原问题的最优值：")
    if h == '1':
        print('z*={}'.format(f(str(-matrix[-1][-1])).limit_denominator()))
    else:
        print('z*={}'.format(f(str(matrix[-1][-1])).limit_denominator()))
def print_duiousulution(matrix):
    print('*' * 20)
    print("对偶问题的最优解：")
    for i in range(n-m,n-1):
        j=i-n+m+1
        print('w{}*={}'.format(j,f(str(-matrix[-1][i])).limit_denominator()),end='，')
def main():
    a = getinput()
    matrix = np.array(a, dtype=np.float64)   #array可以方便地进行整行的操作，而列表可以方便索引
    vect=[]
    h=n
    for i in range(m):
        h=h-i
        vect.append(h)
    pr(matrix,vect)
    while judge(matrix):
        a,matrix,vect = trans(a,matrix,vect)
        pr(matrix,vect)
    print_solution(matrix,vect)
    print_duiousulution(matrix)
if __name__ == '__main__':
	main()

 

4.测试样例