单纯形法与穷举法解线性规划基于python

转载请注明出处

单纯形法

import numpy as np
from numpy import *
from itertools import combinations, permutations
import time

start = time.time()


class Simplex(object):
    def __init__(self, obj, max_mode=False):  # default is solve min LP, if want to solve max lp,should * -1
        self.mat, self.max_mode = np.array([[0] + obj]) * (-1 if max_mode else 1), max_mode

    def add_constraint(self, a, b):
        self.mat = np.vstack([self.mat, [b] + a])

    def _simplex(self, mat, B, m, n):
        while mat[0, 1:].min() < 0:
            col = np.where(mat[0, 1:] < 0)[0][0] + 1  # use Bland's method to avoid degeneracy. use mat[0].argmin() ok?
            row = np.array([mat[i][0] / mat[i][col] if mat[i][col] > 0 else 0x7fffffff for i in
                            range(1, mat.shape[0])]).argmin() + 1  # find the theta index
            if mat[row][col] <= 0: return None  # the theta is ∞, the problem is unbounded
            self._pivot(mat, B, row, col)
        return mat[0][0] * (1 if self.max_mode else -1), {B[i]: mat[i, 0] for i in range(1, m) if B[i] < n}

    def _pivot(self, mat, B, row, col):
        mat[row] /= mat[row][col]
        ids = np.arange(mat.shape[0]) != row
        mat[ids] -= mat[row] * mat[ids, col:col + 1]  # for each i!= row do: mat[i]= mat[i] - mat[row] * mat[i][col]
        B[row] = col

    def solve(self):
        m, n = self.mat.shape  # m - 1 is the number slack variables we should add
        temp, B = np.vstack([np.zeros((1, m - 1)), np.eye(m - 1)]), list(range(n - 1, n + m - 1))  # add diagonal array
        mat = self.mat = np.hstack([self.mat, temp])  # combine them!
        if mat[1:, 0].min() < 0:  # is the initial basic solution feasible?
            row = mat[1:, 0].argmin() + 1  # find the index of min b
            temp, mat[0] = np.copy(mat[0]), 0  # set first row value to zero, and store the previous value
            mat = np.hstack([mat, np.array([1] + [-1] * (m - 1)).reshape((-1, 1))])
            self._pivot(mat, B, row, mat.shape[1] - 1)
            if self._simplex(mat, B, m, n)[0] != 0: return None  # the problem has no answer

            if mat.shape[1] - 1 in B:  # if the x0 in B, we should pivot it.
                self._pivot(mat, B, B.index(mat.shape[1] - 1), np.where(mat[0, 1:] != 0)[0][0] + 1)
            self.mat = np.vstack([temp, mat[1:, :-1]])  # recover the first line
            for i, x in enumerate(B[1:]):
                self.mat[0] -= self.mat[0, x] * self.mat[i + 1]
        return self._simplex(self.mat, B, m, n)


t = Simplex([20, 40, 15, 30, 30, 30, 35, 60, 30, 70])
t.add_constraint([1, 1, 1, 1, 1, 0, 0, 0, 0, 0], 100)
t.add_constraint([1, 0, 0, 0, 0, 1, 0, 0, 0, 0], 20)
t.add_constraint([0, 1, 0, 0, 0, 0, 1, 0, 0, 0], 50)
t.add_constraint([0, 0, 1, 0, 0, 0, 0, 1, 0, 0], 60)
t.add_constraint([0, 0, 0, 1, 0, 0, 0, 0, 1, 0], 30)
t.add_constraint([0, 0, 0, 0, 1, 0, 0, 0, 0, 1], 70)

t.add_constraint([-1, -1, -1, -1, -1, 0, 0, 0, 0, 0], -100)
t.add_constraint([-1, 0, 0, 0, 0, -1, 0, 0, 0, 0], -20)
t.add_constraint([0, -1, 0, 0, 0, 0, -1, 0, 0, 0], -50)
t.add_constraint([0, 0, -1, 0, 0, 0, 0, -1, 0, 0], -60)
t.add_constraint([0, 0, 0, -1, 0, 0, 0, 0, -1, 0], -30)
t.add_constraint([0, 0, 0, 0, -1, 0, 0, 0, 0, -1], -70)
print(t.solve())

end = time.time()
print('以上是单纯形法，运行时间为'+str(end - start)+'s')

穷举法

import numpy as np
from numpy import *
from itertools import combinations, permutations
import time

start = time.time()

C = [20, 40, 15, 30, 30, 30, 35, 60, 30, 70]

A = np.array([[1, 1, 1, 1, 1, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 1, 0, 0, 0, 0],
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0, 1]])
b = np.array([[100], [20], [50], [60], [30], [70]])

lis = list(combinations([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], len(b)))

min = 10000000
c_min = 0
x_min = 0

for i in lis:
    try:
        a = A[:, list(tuple(i))] # 系数矩阵
        c = [C[i] for i in list(tuple(i))] # 目标函数参数
        x = np.linalg.solve(a, b).T[0] # 方程组解
        if (np.array(x) > 0).all():
            if sum([i * j for i, j in zip(x, c)]) < min:
                c_min = list(tuple(i))
                x_min = x
                min = sum([i * j for i, j in zip(x, c)])
    except:
        continue

print(min,c_min,x_min)
end = time.time()
print('以上是穷举法，运行时间为'+str(end - start)+'s')

例子

参考文章

https://blog.csdn.net/lanchunhui/article/details/51824602

https://www.hrwhisper.me/introduction-to-simplex-algorithm/