单纯形法

 

首先需要了解几个概念：

对于一个标准的线性规划问题

设A为m×n满秩的矩阵，将它分块为A=[B N]，其中B为m×m非奇异矩阵，x按B和N的列选择划分成x=（x_B x_N）^T

于是Ax=b可以写成

因此有

取x_N=0，有

这样的解叫做基本解，x_N、x_B分别叫非基本变量和基本变量，N、B分别是基矩阵、非基矩阵

满足非负可行条件的基本解叫做基本可行解

可以证明，每一个基本可行解对应可行域中的一个顶点

 

 

单纯形法是一种广泛应用的线性规划方法，其思想是在顶点上逐个检验最优性，然后寻找相邻的更优顶点，直到确定问题的最优解

在这个算法中，关键的两步是1）检验解是否是最优解，2）寻找下一个基本可行解

这里有几个新的概念：

• 在检验最优解这一步有单纯形乘子、价值系数向量 ，其含义见算法
• 在寻找下一个基本可行解这一步，做法是把原非基矩阵的一列（入基变量）和原矩阵的一列（出基变量）交换

 

单纯形法的步骤如下：

 

下面给出单纯形法的python实现：

# coding=utf-8
from numpy import *


def simple(c, A, b):
    m = min(A.shape)
    n = max(A.shape)
    M = range(0, n)
    E = range(0, m)
    B = A[:, E]
    I = setdiff1d(M, E)
    N = A[:, I]
    i = 1
    x = arange(0, n).astype(float)
    while True:
        if linalg.matrix_rank(B) == m:
            x[E] = dot(linalg.inv(B), b).flatten()
            x[I] = arange(0, n - m)
            if min(x) >= 0:
                break
        E[m - i] += 1
        I = setdiff1d(M, E)
        N = A[:, I]
        B = A[:, E]
        i = i + 1
    # -------
    # E=[2,3,4]
    # I=[0,1]
    # N = A[:, I]
    # B = A[:, E]
    # x=array([0.,0.,3.,2.,16.])
    # -------
    cb = c[E]
    cn = c[I]
    # 单纯形乘子
    y = dot(linalg.inv(B.T), cb)
    # 简约价值系数
    cv = cn - dot(N.T, y)
    eps = 10 ** (-7)
    r = b
    while min(cv) < 0:
        p = I[where(cv == min(cv))[0][0]]
        ap = dot(linalg.inv(B), A[:, p])
        bq = r[ap > eps] / ap[ap > eps]
        bq = min(bq)
        q = E[where(r / ap == bq)[0][0]]
        out = where(array(E) == q)[0][0]
        enter = where(array(I) == p)[0][0]
        E[out] = p
        I[enter] = q
        N = A[:, I]
        B = A[:, E]
        x[E] = dot(linalg.inv(B), b).flatten()
        r = x[E]
        x[I] = zeros(n - m)
        cb = c[E]
        cn = c[I]
        y = dot(linalg.inv(B.T), cb)
        cv = cn - dot(N.T, y)
    print(x)


c = array([-2, -3, 0, 0, 0])
A = array([[-1, 1, 1, 0, 0], [-2, 1, 0, 1, 0], [4, 1, 0, 0, 1]])
b = array([3, 2, 16])
simple(c, A, b)

运行结果：

即x₁=2.6 x₂=5.6 x₃=0 x₄=1.6 x₅=1 是该问题的最优解