【运筹学】单纯形法（两阶段法）Python实现（第二版）

1 单纯形法（两阶段法）--Python

        两阶段法主要是用于初始基可行解难以寻找的情况；

1.1 创建第一阶段单纯形表，添加人工变量

import numpy as np
import sys
def create_first_tableau(A,b,num_constraints,num_variables):
    basic_vars = np.arange(num_variables, num_variables + num_constraints)
    non_basic_vars = np.arange(num_variables)
    # np.hstack水平组合向量或矩阵；array.reshape(-1,1)将b转化为n行1列的矩阵；
    tableau=np.hstack((A, b.reshape(-1,1) ))
    for i in range(num_constraints): # 保持b列始终非负
        if tableau[i,-1]<0:
            tableau[i,:]*=-1
    basic_matrix=np.eye(num_constraints)
    tableau=np.insert(tableau,num_variables,basic_matrix,axis=1)
    c_row=np.zeros(num_constraints+num_variables+1)# 构建系数行
    for i in range(num_constraints):
        c_row-=tableau[i,:]
    c_row[basic_vars]=0
    tableau=np.vstack((tableau,c_row))
    print("第一阶段：\n")
    print("Iteration 0：\n",tableau)
    print("Basic Variables:", basic_vars+1)
    print("Non-Basic Variables:", non_basic_vars+1)
    return tableau,basic_vars,non_basic_vars

1.2 创建第二阶段单纯形表，消除人工变量

def create_second_tableau(tableau,c,num_constraints,num_variables,basic_vars,non_basic_vars):
    artificial_vari=np.arange(num_variables, num_variables + num_constraints)
    non_basic_vars=non_basic_vars[~np.isin(non_basic_vars,artificial_vari)]
    c_row=np.hstack((c,np.zeros(num_constraints+1))) # 构建系数行
    tableau[-1,:]=c_row
    new_tableau=[i for i in range(num_variables)]
    new_tableau.append(-1) # 消除人工变量列
    tableau=tableau[:,new_tableau]
    for i,arg in enumerate(basic_vars):
        tableau[-1,:]+=-tableau[-1,arg]*tableau[i,:]
    print('\n第二阶段：\n')
    print("Iteration 0：\n",tableau)
    print("Basic Variables:", basic_vars+1)
    print("Non-Basic Variables:", non_basic_vars+1)
    return tableau,basic_vars,non_basic_vars

1.3 定义换基(迭代)过程

# 首先判断目标函数系数向量中非基变量的值是否非负，若全是非负数，则达到最优；若存在负数，则取最小值对应的变量作为入基变量；
# 原理是利用目标函数系数向量*变化量向量（改进方向/非基变量列），观察是否减小；
# 由于目标函数系数向量跟着变化（单纯形字典），用非基变量表示目标函数，而换基时是固定一个非基变量，变化一个非基变量，因此可以只看目标函数系数即可
def pivot(tableau,basic_vars,non_basic_vars):
    tableau_N=tableau[-1,non_basic_vars]
    arg_min_N=np.argmin(tableau_N)
    entering_vari_col=non_basic_vars[arg_min_N]
    if tableau[-1,entering_vari_col]>=0:
        print('当前解为最优解') 
        return False
    else:
        print(f'入基变量为x{entering_vari_col+1};',end='')
        step=np.where(tableau[:-1,entering_vari_col]>0)[0]
        if len(step)==0:
            print('\n此LP问题无界')
            sys.exit()
        else:
            #最小比值法确定步长和出基变量；存在疑问：基变量的变化量会为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) #需要取正值最小值的索引
            print(f'出基变量为x{basic_vars[leaving_vari_row]+1};')
            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)
            convert=basic_vars[leaving_vari_row]
            basic_vars[leaving_vari_row]=non_basic_vars[arg_min_N]
            non_basic_vars[arg_min_N]=convert
            print("Basic Variables:", basic_vars+1)
            print("Non-Basic Variables:", non_basic_vars+1)
            return True

1.4 提取目标函数值和解的情况

def extract_optimal(tableau,basic_vars,non_basic_vars):
    optimal_value = -tableau[-1, -1]  # tableau[-1, -1]是目标函数最小值的负数
    optimal_x=np.zeros(len(tableau[0])-1)
    for i,arg in enumerate(non_basic_vars):
        optimal_x[arg]=0
    for i,arg in enumerate(basic_vars):
        optimal_x[arg]=tableau[i,-1]
    print("目标函数最小值:", optimal_value)
    print("解的情况：",end='')
    for i,val in enumerate(optimal_x):
        print(f'x{i+1}={val}',end='; ')
    return optimal_x,optimal_value

2 输入参数等

# 输入的参数一定是标准化后的线性规划问题的参数
c = np.array([-2,-1,-9,0,0])  # 目标函数的系数
A = np.array([[1,1,0,1,0],
              [-2,0,1,0,0],
              [0,3,5,0,-1]])  # 约束矩阵
b = np.array([18, -2, 15])  # 约束的右端项
num_constraints,num_variables=A.shape

2.1 求解第一阶段

tableau,basic_vars,non_basic_vars=create_first_tableau(A,b,num_constraints,num_variables)
iteration=1
while True:
    print(f"\nIteration {iteration}:",end='')
    iteration+=1
    if pivot(tableau,basic_vars,non_basic_vars)==False:
        break
    elif iteration==100: # 防止进入循环
        print('\n此问题无可行解')
        sys.exit()
optimal_x1,optimal_value1=extract_optimal(tableau,basic_vars,non_basic_vars)
if abs(optimal_value1)<=0.000001: # 该例题中目标函数值趋近于0，不知道问题所在
    print("\n人工变量之和为0，原问题存在可行解，进行第二阶段计算")
else:
    print("\n人工变量之和不为0，原问题不存在可行解，停止计算")
    sys.exit()

2.2 求解第二阶段

tableau,basic_vars,non_basic_vars=create_second_tableau(tableau,c,num_constraints,num_variables,basic_vars,non_basic_vars)
iteration=1
while True:
    print(f"\nIteration {iteration}:",end='')
    iteration+=1
    if pivot(tableau,basic_vars,non_basic_vars)==False:
        break
    elif iteration==100:
        print('\n此问题无可行解')
        sys.exit()
optimal_x2,optimal_value2=extract_optimal(tableau,basic_vars,non_basic_vars)

3 例题及运行结果

3.1 例题及标准化过程

3.2 算法运行结果

3.3 验证结果（Lingo）

4 修改说明 

        昨天的代码中没有注意到消除人工变量后非基变量索引的变化，作此修改non_basic_vars=non_basic_vars[~np.isin(non_basic_vars,artificial_vari)]，并规范了更多代码运行流程；仍然存在的问题：1.没有将线性规划问题标准化的过程，需要读者先将问题进行标准化；2.仍未考虑退化解以及循环的情况；

5 关于作者

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