matlab单纯形法

%单纯形法
function [xb,z] = text3(a,b0,c)
%赋初值

[m,n] = size(a);
f = (c==0);
number = sum(f(:));
e = eye(number);
B = e;
x = [1:n;ones(1,n)];
xb = zeros(m,1);%基变量标记
xn = [];%非基变量
cn = [];%非基变量的价格系数
cb = zeros(1,m);
xin = 0;xout = 0;%换入变量和换出变量

%找到单位矩阵
%处理基变量以及对应数值
for j = 1:m
    for i =1:n
        if a(:,i) == e(:,j)
            xb(j) = i;
            x(2,i) = 0; %将x中1变为0，变为基变量
            cb(j) = c(i);%基变量的价值系数
        end
    end
end
%处理非基变量以及对应数值
for i = 1:n
    if x(2,i) == 1
        xn = [xn;x(1,i)];%得到非基变量
        cn = [cn,c(i)];%得到非基变量的价格系数
    end
end

b = b0;
jianyan = c;%检验数的初值设为价格系数向量c
jianyanmax = max(jianyan);
kk = 1;

while(jianyanmax>0.000000000001)
    in = find(jianyan == jianyanmax);
    in = in(1);
    xin = in;
    P = B*a(:,in);
    xita = 1000*ones(number,1);%设定θ初值
    
    for i = 1:m
        if sum(P<=0)==number && jianyan(i)>0
            disp('LP问题为无界解')
            xb = [];
            z = [];
            return
        end
        if P(i)>0%判断θ是否符合
            xita(i) = b(i)/P(i);%改变θ的值，由此确定换出变量
        end
    end
    
    [~,out] = min(xita);
    out = out(1);
    xout = xb(out);%得到换出变量
    xb(out) = xin;%将换入变量换入基变量中
    cb(out) = c(xin);%cb为基变量的价格系数
    weizhi = find(xn==xin);
    xn(weizhi) = xout;%xn为非基变量的位置
    cn(weizhi) = c(xout);%cn为非基变量的价格系数
    %单纯形法的矩阵表示：
    EE = e;
    EE(:,out) = P;
    E = inv(EE);
    B = E*B;
    bb = b0;
    b = B*bb;
    jianyan = c-cb*B*a;
    jianyanmax = max(jianyan);
    kk = kk+1;
end
z = cb * b;
xb = [xb b];
end