Python解决两阶段法规划问题

前言

在做运筹学作业的过程中，我被复杂的数学逻辑困扰。于是有了做一个程序帮助我写作业的想法。
亮点：借助了前人的成果，改进的地方在于展示了每一步迭代的详细过程，输出美观的表格，贴近实际解题。改进了输入的逻辑,傻瓜式操作，仅需按提示输入条件。觉得有帮助可以点个赞。

代码

import numpy as np
import time

class Simplex(object):
    def __init__(self, c, A_ub=None, b_ub=None, A_eq=None, b_eq=None):
        self.c = c  # 系数矩阵
        self.A_ub = A_ub if A_ub is not None else np.array([])
        self.b_ub = b_ub if b_ub is not None else np.array([])
        self.A_eq = A_eq if A_eq is not None else np.array([])
        self.b_eq = b_eq if b_eq is not None else np.array([])
        self.N_x = 0  # 变量个数
        self.meq = 0  # 等式约束个数
        self.mub = 0  # 不等式约束个数
        self.m = 0  # 约束总数
        self.T = []  # 单纯形表
        self.sign = []  # 基变量索引
        self.F = 1  # 迭代标志（1: 继续，0: 无界）

    def Initial_value(self):
        self.N_x = len(self.c)
        # 确定约束个数
        if self.b_eq.any():
            self.meq = len(self.b_eq)
        if self.b_ub.any():
            self.mub = len(self.b_ub)
        self.m = self.mub + self.meq  # 约束总数
        # 创建初始单纯形表
        self.T = np.zeros([self.m + 1, self.N_x + self.m + 1], dtype='float')
        b = self.T[:-1, -1]
        if self.meq > 0:
            self.T[:self.meq, :self.N_x] = self.A_eq
            b[:self.meq] = self.b_eq
        if self.mub > 0:
            self.T[self.meq:self.m, :self.N_x] = self.A_ub
            b[self.meq:self.m] = self.b_ub
        np.fill_diagonal(self.T[:self.m, self.N_x:self.N_x + self.m], 1)
        self.T[-1, self.N_x:self.N_x + self.meq] = -1 * np.ones(self.meq)
        for i in range(0, self.meq):
            self.T[-1] += self.T[i]
        self.sign = list(range(self.N_x, self.N_x + self.m))

    def solve(self):
        num = 0
        flag = True
        while flag:
            print(f"\n第 {num} 次迭代的单纯形表：")
            self.print_table()
            if max(list(self.T[-1][:-1])) <= 0:
                flag = False
            else:
                num += 1
                self.F = self.calculate()
            if self.F == 0:
                break

    def calculate(self):
        H = list(self.T[-1, :-1])
        j_num = H.index(max(H))
        D = []
        for i in range(0, self.m):
            if self.T[i][j_num] == 0:
                D.append(float("inf"))
            else:
                D.append(self.T[i][-1] / self.T[i][j_num])
        if max(D) <= 0:
            return 0
        i_num = D.index(min([x for x in D if x >= 0]))
        self.sign[i_num] = j_num
        t = self.T[i_num][j_num]
        self.T[i_num] /= t
        for i in [x for x in range(0, self.m + 1) if x != i_num]:
            self.T[i] -= self.T[i_num] * self.T[i][j_num]
        return 1

    def change(self):
        self.T[:, self.N_x:self.N_x + self.meq] = 0
        self.T[-1, 0:self.N_x] = -self.c
        for i in range(0, self.m):
            self.T[-1] -= self.T[i] * self.T[-1][int(self.sign[i])]

    def Main(self):
        self.Initial_value()
        print("初始单纯形表：")
        self.print_table()
        if self.meq > 0:
            print("第一阶段")
            self.solve()
            self.change()
            print("第二阶段")
            self.solve()
        else:
            print("单阶段")
            self.change()
            self.solve()
        if self.F == 1:
            print("最优解：")
            j = 0
            for i in range(0, self.N_x):
                if i not in self.sign:
                    print(f"x{i + 1} = 0")
                else:
                    print(f"x{i + 1} = {self.T[j][-1]}")
                    j += 1
            print("最优值：\n", self.T[-1][-1])
        else:
            print("出现无界解")

    def print_table(self):
        header = ["变量"] + [f"x{i+1}" for i in range(self.N_x + self.m)] + ["常数"]
        rows = []
        for i, row in enumerate(self.T):
            if i < self.m:
                rows.append([f"x{self.sign[i]+1}"] + list(np.round(row, 2)))
            else:
                rows.append(["Z"] + list(np.round(row, 2)))
        print("\t".join(header))
        for row in rows:
            print("\t".join(map(str, row)))

def user_input():
    print("请输入目标函数的系数 (例如：6x1 + 3x2 + 4x3，输入 6, 3, 4)：")
    c = list(map(float, input().strip().split(',')))

    print("请输入等式约束条件 (例如：x1 + x2 + x3 = 120，输入 1, 1, 1 = 120)：")
    A_eq = []
    b_eq = []
    while True:
        constraint = input("等式约束 (输入 'done' 完成输入)：")
        if constraint.lower() == 'done':
            break
        lhs, rhs = constraint.split('=')
        A_eq.append(list(map(float, lhs.strip().split(','))))
        b_eq.append(float(rhs.strip()))

    print("请输入不等式约束条件 (例如：x1 >= 30 转换为 -1x1 <= -30，输入 -1, 0, 0 <= -30)：")
    A_ub = []
    b_ub = []
    while True:
        constraint = input("不等式约束 (输入 'done' 完成输入)：")
        if constraint.lower() == 'done':
            break
        lhs, rhs = constraint.split('<=')
        A_ub.append(list(map(float, lhs.strip().split(','))))
        b_ub.append(float(rhs.strip()))

    return np.array(c), np.array(A_ub), np.array(b_ub), np.array(A_eq), np.array(b_eq)

if __name__ == '__main__':
    c, A_ub, b_ub, A_eq, b_eq = user_input()

    time_start = time.perf_counter()
    problem = Simplex(c, A_ub, b_ub, A_eq, b_eq)
    problem.Main()
    time_end = time.perf_counter()
    print('运行时间: %.2f 秒' % (time_end - time_start))