[Python知识库]python实现单纯形法迭代形式（待更新人工法+对偶）

输入格式与matlab的linprog公式类似链接，之前还写过这个已经忘了，可以参考，但是这个的写法要求不等式必须是≤，有需求可以自己变通

单纯形法的迭代公式建议看厦大运筹学之规划论

剩下的实现和原理之后再写，这个确实挺好用的🐶
我把例子找了三个朴素的例子做了实现，结果还可以

import numpy as np


def LinearProgram(f, A, b, Aeq=-1, beq=-1):
    """
    传入matlab样式参数，通过判断选择单纯形法，选择性添加人工法，对偶问题求解
    :param f:   目标函数
    :param A:   不等式约束系数矩阵
    :param b:   不等式约束常数矩阵
    :param Aeq: 等式约束系数矩阵
    :param beq: 等式约束常数矩阵
    :return: 若有值则返回函数值与对应基变量
    """

    # 只有不等式约束采用单纯形法基础方法解决，先进行
    if Aeq == -1:
        # 构造基向量组
        basic_vector = np.eye(A.shape[0])
        sheet = np.hstack((A, basic_vector))
        # 标记
        basic_index = [i for i in range(A.shape[1], A.shape[0] + A.shape[1])]
        # 连接
        f = np.append(f, np.zeros((1, A.shape[0])))
        ans = SimplexMethod(f, sheet, b, basic_index)
        if ans[0] == 3:
            print(ans[1])
        else:
            x_variable = np.zeros((f.shape[0], 1))
            for i in range(len(ans[2])):
                x_variable[ans[2][i]] = ans[3][i]
            fval = np.dot(f, x_variable)
            if ans[0] == 1:
                print("{0};\n函数取值为{1};\n最终函数值为{2}.\n"
                      .format(ans[1], ' '.join(map(str, x_variable)), fval))
            else:
                print("{0};\n函数取值为{1};\n最终函数值为{2}.\n"
                      .format(ans[1], ' '.join(map(str, x_variable)), fval))

    return 1, 2


def SimplexMethod(f, sheet, b, index):
    """
    朴素单纯形法求解
    :param f: 处理过的函数项
    :param sheet: 表
    :param b: 常数项
    :param index: 基解变量index
    :return: 结果
    """
    basic_c = f[index]
    while True:
        # 生成zeta
        zeta = np.full(sheet.shape[1], np.nan)
        for i in range(f.shape[0]):
            if i not in index:
                zeta[i] = f[i] - np.dot(basic_c, sheet[:, i])
        # 找出进基下标
        push_in_index = np.nanargmin(zeta)

        if zeta[push_in_index] > 0:
            return [1, "唯一最优解", index, b]
        elif zeta[push_in_index] >= 0:
            return [2, "无穷最优解", index, b]

        # 生成theta
        theta = np.full(len(index), np.nan)
        for i in range(len(index)):
            if sheet[i, push_in_index] > 0:
                theta[i] = b[i] / sheet[i, push_in_index]
        # 找出出基下标
        try:
            push_out_index = np.nanargmin(theta)
        except ValueError:
            return [3, "无最优解/有无界解"]

        # 换基
        index[push_out_index] = push_in_index
        basic_c[push_out_index] = f[push_in_index]

        # 自身单位化
        unit_change = sheet[push_out_index, push_in_index]
        sheet[push_out_index] /= unit_change
        b[push_out_index] /= unit_change

        # 整体单位化
        for i in range(sheet.shape[0]):
            if i == push_out_index:
                continue
            full_change = sheet[i, push_in_index]
            sheet[i] -= full_change * sheet[push_out_index]
            b[i] -= full_change * b[push_out_index]


if __name__ == '__main__':
    f = np.array([-1, -2])
    A = np.array([[1, 0], [0, 1], [1, 2]])
    b = np.array([4, 3, 8])
    # Aeq = np.array([])
    # beq = np.array([])
    LinearProgram(f, A, b)

    f = np.array([-6, -4])
    A = np.array([[2, 3], [4, 2]])
    b = np.array([100, 120])
    # Aeq = np.array([])
    # beq = np.array([])
    LinearProgram(f, A, b)

    f = np.array([-2, -1])
    A = np.array([[-1, 1], [2, -5]])
    b = np.array([5, 10])
    # Aeq = np.array([])
    # beq = np.array([])
    LinearProgram(f, A, b)