simplex

import numpy as np
#首先给出标准化矩阵A及右端b，目标函数系数矩阵c
def baseVarShow(base_var,b,c):
    var_list = ""
    z = 0
    b = b.flatten()
    for i in range(len(b)):
        var_list += f'X{base_var[i] + 1} = {b[i]}；'
        z += c[int(base_var[i])] * b[i]
    print("基本可行解变量分别为(未显示的为非基变量，令其为零)：")
    print(var_list)
    print("目标函数值为：")
    print(z)
    pass
def findOutBase(base_var, p, b):
    num = np.inf
    out_base = 0
    for i in range(len(p)):
        if p[i] > 0:
            if num >= b[i]/p[i] and b[i]/p[i] >= 0:
                num = b[i]/p[i]
                out_base = base_var[i]
    return out_base
def  algorithm(A, b, c, base_var, F):
    c_base = c[base_var]#获取基本可行解的系数
    mat_A = np.mat(A)
    mat_B = np.mat(A[:, [base_var]])#获取基矩阵
    sigma_list = c - c_base * mat_B.I * mat_A
    if sum(sigma_list.T <= 0) == np.shape(sigma_list)[1]:#判断是否为最优解
        print("#####################获取最优解：#####################")
        baseVarShow(base_var, b, c)
        np.savetxt("D:/作业/第一阶段结束matrix_A.txt", mat_A, fmt = '%.5f')
        np.savetxt("D:/作业/第一阶段结束vector_b.txt", b, fmt = '%.5f')
        np.savetxt("D:/作业/第一阶段结束base_var.txt", base_var, fmt = '%.0f')
        print("此时矩阵A为：")
        print(A.astype(float))
    else:#若不是最优解则进入递归
        print("")
        if F == '1':
            choose = input("是否查看当前迭代状态，是则输入1：")
            if choose == '1':
                print("################################当前状态为：################################\n")
                baseVarShow(base_var, b, c)
                print("此时矩阵A为：")
                print(A.astype(float))
                print("目标系数矩阵为：", c)
                print("检验数矩阵为：", sigma_list.astype(float))
                print("")
                print("############################################################################")
        in_base = np.argwhere(sigma_list == max(sigma_list.T))[0][1]#选择进基变量·
        p = A[:, in_base]
        if sum(p <= 0) == len(p):
            print(f"当选择X{in_base+1}做进基变量时，出现对应向量各数值均小于等于0，无法利用最小比值准则，故无最优解")
            print("此时对应基向量为：", p)
            return
        out_base = findOutBase(base_var, p, b)#获取离基变量
        base_var[base_var.index(out_base)] = in_base
        mat_B1 = np.mat(A[:, [base_var]]).I
        algorithm(np.array(mat_B1 * mat_A), np.array(mat_B1 *  np.reshape(b,  (len(b), 1))), c, base_var, F)

def simplexMethod(A, b, c, F, base_var = []):#模拟单纯型表算法
    row = []
    flag = 0
    if len(base_var) == 0:
        flag = 1
    rank_A = np.linalg.matrix_rank(A)
    if rank_A == np.shape(A)[0]:#若矩阵A的秩不等于A的行数则说明矩阵输入有误
        for i in range(np.shape(A)[1]):
            if sum(A[:, i] == 1) == 1:
                if sum(A[:, i] == 0) == rank_A - 1 and int(np.argwhere(A[:, i] == 1)) not in row:
                    if flag:
                        base_var.append(i)#获取初始基本可行解
                        row.append(int(np.argwhere(A[:, i] == 1)))
        if len(base_var) == rank_A:
            print("#####################初始基本可行解：#####################")
            baseVarShow(base_var, b, c)
            algorithm(A, b, c, base_var, F)
        else:
            print("请检查输入矩阵是否为标准化矩阵")
    else:
        print("输入有误")
def bigM(A,b,c):
    print("======================================大M法计算开始：======================================")
    M = 999
    print("M设为999")
    n1 = np.shape(A)[1]
    n2 = np.shape(A)[0]
    x = ""
    for i in range(n2):
        x += f"X{i+n1+1}；"
    print("新增人工变量：", x)
    eyes = np.eye(np.shape(A)[0])
    A = np.append(A, eyes, axis = 1)
    c = np.append(c, -np.ones(np.shape(A)[0])*M, axis = 0)
    print("")
    F = input("是否查看迭代过程，查看请输入1：")
    simplexMethod(A,b,c,F)
    print("======================================大M法计算结束：======================================")
def twoStep(A,b,c):
    print("======================================二阶段法第一阶段开始：======================================")
    n1 = np.shape(A)[1]
    n2 = np.shape(A)[0]
    x = ""
    for i in range(n2):
        x += f"X{i+n1+1}；"
    print("新增人工变量：", x)
    eyes = np.eye(np.shape(A)[0])
    c1 = np.append(np.zeros(np.shape(A)[1]), -1*np.ones(np.shape(eyes)[0]), axis = 0)
    A = np.append(A, eyes, axis = 1)
    print("")
    F = input("是否查看迭代过程，查看请输入1：")
    simplexMethod(A,b,c1,F)
    print("======================================二阶段法第一阶段结束：======================================")
    print("")
    print("======================================二阶段法第二阶段开始：======================================")
    A1 = np.loadtxt("D:/作业/第一阶段结束matrix_A.txt")
    b1 = np.loadtxt("D:/作业/第一阶段结束vector_b.txt")
    base_var = np.loadtxt("D:/作业/第一阶段结束base_var.txt").astype(int).tolist()
    print("")
    F = input("是否查看迭代过程，查看请输入1：")
    simplexMethod(A1[:,:n1], b1, c, F, base_var)
    print("======================================二阶段法第二阶段结束：======================================")

print("##################################################################################################")
print("使用说明：\n请建立以下文件夹及相关文件：\nD:/作业/初始标准化矩阵matrix_A.txt\nD:/作业/初始标准化矩阵vector_b.txt\nD:/作业/初始标准化矩阵vector_c.txt\n"
      "并分别在matrix_A中输入标准形式约束矩阵A，在vector_b中输入标准形式约束右端向量b，在vector_c中输入目标系数矩阵c")
print("##################################################################################################")
choose = input("输入1选择大M，否则二阶段：")
if choose == '1':
    A = np.loadtxt("D:/作业/初始标准化矩阵matrix_A.txt")
    b = np.loadtxt("D:/作业/初始标准化矩阵vector_b.txt")
    c = np.loadtxt("D:/作业/初始标准化矩阵vector_c.txt")
    bigM(A,b,c)
else:
    A = np.loadtxt("D:/作业/初始标准化矩阵matrix_A.txt")
    b = np.loadtxt("D:/作业/初始标准化矩阵vector_b.txt")
    c = np.loadtxt("D:/作业/初始标准化矩阵vector_c.txt")
    twoStep(A,b,c)
input()