Python单纯形法求解线性规划 用单纯形法解线性规划

单纯形法是由George Dantzig于1947年提出的一种求解线性规划问题的算法。在线性规划问题中，如果最优解存在，那么它一定可以在可行区域的顶点中找到。单纯形法的基本思路是：先找出可行域的一个顶点，根据一定规则判断其是否最优；若否，则转换到与之相邻的另一顶点，并使目标函数值更优；如此循环，直到找到最优解为止。

在单纯形法中，凸集的概念尤为重要。凸集上求最优解是凸优化的一个分支，对于简化问题具有重要意义。而单纯形可以理解为一个凸集，在N维空间中，它是由N+1个顶点的凸包组成的多胞体。

单纯形法的具体步骤如下：

1. 求出线性规划的初始基可行解，并列出初始单纯形表。这通常涉及到基于约束条件方程组的系数矩阵，通过寻找或构造单位矩阵的方法，确定基变量，从而求出初始基本可行解。再利用初始基本可行解及线性规划模型提供的信息，编制初始单纯形表。
2. 进行最优性检验。将检验数cj-zj作为判断基本可行解是否为最优解的标准。若所有非基变量的检验数cj-zj<0，则已经达到最优解，计算停止。
3. 从一个基可行解转换到另一个目标函数值更大的基可行解，列出新的单纯形表。这涉及到确定换入变量和换出变量，并进行相应的置换操作。换入变量通常选择最大检验数对应的非基变量，而换出变量则通过一定的计算规则确定。完成换入和换出后，需要更新单纯形表，包括变更目标函数系数和进行行运算等。
4. 重复第二、三步，直到计算终止。这意味着要么找到了最优解，要么判断出线性规划问题没有最优解。

在整个过程中，单纯形法通过不断迭代，从可行域中的一个点开始，逐步找到目标函数值更优的顶点，直到找到最优解或确定问题无解。

需要注意的是，单纯形法是一种比较复杂的算法，需要一定的数学基础和计算能力。在实际应用中，通常会使用计算机程序来辅助完成计算过程。下面就用MATLAB来实现单纯型法：

c=[-2 1 -3 1 0 0];
A=[1 3 1 1 0 0;
   -2 1 1 0 1 0;
   -1 6 -1 0 0 1;];
b=[6;3;4];
[ x,cx,B,k ] = simplex( c,A,b )

function [x,cx,B,k] = simplex (c,A,b)
% 线性规划典式对应的单纯形法
% c是目标函数中x的系数行向量，A是约束条件Ax=b中x的系数矩阵,b是右端项
% 本程序要求A由（N,E）构成，其中E是单位阵，b>=0,即x中的后m个分量是基变量，此类问题易于找到可行基E
% LP问题要求变量x非负
% 输出最优解x，最优值cx，最优基B，换基迭代次数k

%构建初始单纯形表AA
[m,n]=size(A);
xb=n-m+1:n;
AA=zeros(m+1,n+1);
for j=1:n-m
    AA(1,j+1)=c(j)-c(xb)*A(:,j);
end
AA(1,1)=-c(xb)*b;
AA(2:m+1,1)=b;
AA(2:m+1,2:n+1)=A;

%B=zeros(m,m);
%B=A(:,xb);
k=0;
x0=zeros(n,1);
while 1
    % 当前基本可行解
x0(xb)=AA(2:m+1,1);

%选择进基变量
q=0;
for j=1:n
    if AA(1,j+1)<0
        q=j;
        break;
    end
end
if q==0    %满足最优性条件
    x=x0;
    B=A(:,xb);
    cx=-AA(1,1);
    return;
else
    %选择离基变量
    theta=+Inf;
    for i=1:m
        if AA(i+1,q+1)>0 & AA(i+1,1)/AA(i+1,q+1)<theta
            theta=AA(i+1,1)/AA(i+1,q+1);
            p=i;
        end
    end
    if theta==Inf
        disp('This LP does not have a finite solution!');
        x=Inf;
        cx=-Inf;
        B=[];
        return;
    else
        xb(p)=q;
        %进行换基运算，初等行变换
        AA(p+1,:)=AA(p+1,:)/AA(p+1,q+1);
        for i=1:m+1
            if i~=p+1
                aa=AA(i,q+1);%注意此处先要保存元素AA(i,q+1)，否则该元素将在换基运算中变成0，导致在其后无法执行相应运算
                for j=1:n+1
                    AA(i,j)=AA(i,j)-AA(p+1,j)*aa;
                end
            end
        end
    end
   
end
%AA %可以获取每一步换基运算后的单纯形表
k=k+1;
%B=A(:,xb); %可以获取每一步换基后的基矩阵
end

end

运行结果如下图，可以看出上述代码是十分准确的：