【运筹学】学习笔记--单纯形法Python实现

1 利用Python单纯形法

1.1 创建单纯形表

import numpy as np
# 本例中需要采用标准型LP问题（添加人工变量、目标函数最小值、变量取值非负）
# 单纯形表前m行n列为系数表，m+1行为目标函数系数，n+1行为约束值，-(m+1,n+1)为目标函数值；
# 注：此案例中默认添加的人工变量系数矩阵为初始基矩阵；目前没有出基变量的索引（相当于没有解的索引）；
def crate_tableau(A,b,c):
    # np.hstack水平组合向量或矩阵；array.reshape(-1,1)将b转化为n行1列的矩阵；
    tableau=np.hstack((A, b.reshape(-1,1) ))
    # 需要增加判断，是求最小值还是最大值；
    c_row=np.hstack((c, np.zeros(1)))
    tableau=np.vstack((tableau,c_row))
    print('Iteration 0：\n',tableau)
    return tableau

1.2 定义换基过程

# 首先判断目标函数系数向量中非基变量的值是否非负，若全是非负数，则达到最优；若存在负数，则取最小值对应的变量作为入基变量；
# 原理是利用目标函数系数向量*变化量向量（改进方向/非基变量列），观察是否减小；
# 由于目标函数系数向量跟着变化（单纯形字典），用非基变量表示目标函数，而换基时是固定一个非基变量，变化一个非基变量，因此可以只看目标函数系数即可
def pivot(tableau):
    entering_vari_col=np.argmin(tableau[-1,:-1]) #取系数行中最小的值
    if (tableau[-1,entering_vari_col]>=0):
        print('当前解为最优解') 
        return False
    else:
        print(f'入基变量为x{entering_vari_col+1};',end='')
        #最小比值法确定步长和出基变量；存在疑问：基变量的变化量会为0吗？如果变化为0，则约束没变
        ratios=tableau[:-1,-1]/tableau[:-1,entering_vari_col] # b列是基变量的取值，始终非负
        ratios=np.where(ratios<0,np.inf,ratios)
        leaving_vari_row=np.argmin(ratios) #需要取非负最小值的索引
        if ratios[leaving_vari_row]<0:
            print('此LP问题无界')
            return False
        else:
            print(f'出基变量为x{leaving_vari_row+3};')
            tableau[leaving_vari_row,:]/=tableau[leaving_vari_row,entering_vari_col]
            for i in range(len(tableau)):
                if i!=leaving_vari_row:
                    tableau[i,:]-=tableau[i,entering_vari_col]*tableau[leaving_vari_row,:]
            print(tableau)
            return True

1.3 提取最优解和解的情况

def optimal(tableau,num_variable):
    optimal_value = -tableau[-1, -1]  # tableau[-1, -1]是目标函数最小值的负数
    optimal_x=np.zeros(num_variable)
    for i in range(len(A)):
        array=np.where(tableau[i,:num_variable]==1)
        if len(array) == 1:
            optimal_x[array[0]]=tableau[i,-1]
    return optimal_value,optimal_x

2 输入参数等

num_variable=2
c = np.array([-12, -9, 0, 0, 0, 0])  # 目标函数的系数
A = np.array([[1, 0, 1, 0, 0, 0],
              [0, 1, 0, 1, 0, 0],
              [1, 1, 0, 0, 1, 0],
              [4, 2, 0, 0, 0, 1]])  # 约束矩阵
b = np.array([10, 15, 17.5, 48])  # 约束的右端项
tableau=crate_tableau(A,b,c)
iteration=1
while True:
    print(f"\nIteration {iteration}:",end='')
    iteration+=1
    if pivot(tableau)==False:
        break

opt_value,solution=optimal(tableau,num_variable)
print("\n最优解：")
print("目标函数最小值:", opt_value)
for i,val in enumerate(solution):
    print(f'x{i}={val}',end='; ')

 3 问题实例及运行结果

4 关于问题 

        本文旨在记录作者学习运筹学过程中复现经典算法的一步步过程，可能存在算法流程不完善、表述不清楚的情况，请不吝您的宝贵意见，作者会持续学习和完善！