originalsimpleM.m

 

function [x,fval,flag,iteration]=originalsimpleM(C,A,b)

 

%原始单纯形法（大M法，无需给出初始基变量）

%Programmed by Liyang(faruto's Studio~!) BNU MATH

%Email:liyangbnu@mail.bnu.edu.cn  QQ:516667408

%last modified 2008.4.27

%求解标准型线性规划:max C*X;s.t. A*X=b (b>=0);X>=0;

%输入：C是n维行向量，A是m*n的系数矩阵，b是m维列向量

%输出：x最优解（如果有的话），fval最优值，flag解的状态说明，interation求解时的循环次数

%flag 最终解的状态说明：

%flag = 0   LP converged to a solution x

%flag = 1   Inf feasible solutions

%flag = 2   LP is unbounded

%flag = 3   No feasible point was found

 

M = 1e5;

[m,n] = size(A);

A = [A,eye(m)];

for run = 1:1:m

    C = [C,-M];

end

 

XB = (n+1:n+m);%XB承装初始基变量的下标

n = n+m;

 

iteration = 0;

while 1

    iteration = iteration +1;

    flag = 0;

   

%%

    sigma = zeros(1,n);

    for col = 1:n

        temp = 0;

        for row = 1:m

            temp = temp + C(XB(row))*A(row,col);

        end

        sigma(col) = C(col)-temp;

    end

   

%%

    if sigma<=0  

       

        x = zeros(1,n);

        for row = 1:m

            x(XB(row)) = b(row);

        end

        fval = C*x';

       

        for row = 1:m

            for temp = n-m+1:n

                if XB(row)==temp && x(temp)~=0

                    flag = 3;

                    x = 0;

                    fval = 0;

                    break;

                end

            end

            if flag == 3

                break;

            end

        end

        if flag == 3

            break;

        end

       

        for col = 1:n

            tflag = 0;

            for row = 1:m

                if col == XB(row)

                    tflag = 1;

                    break;

                end

            end

            if tflag == 0

                if sigma(col) == 0

                    flag = 1;                  

                    break;

                end

            end

        end

        if flag == 1

            x = x(:,1:n-m);

            break;

        end

       

        if flag == 0; 

            x = x(:,1:n-m);

            fval;

            break;

        end

       

    else

        for col = 1:n

            if sigma(col)>0 & A(:,col)<=0

                flag = 2;

                x = 0;

                fval = 0;

                break;

            end

        end

        if flag == 2

            break;

        end

       

%%

        temp = 0;

        for col = 1:n

            if sigma(col)>temp

                temp = sigma(col);

                intobase = col;%入基变量的下标

            end

        end

       

        theta = zeros(1,m);

        for row = 1:m

            if A(row,intobase)>0

                theta(row) = b(row)/A(row,intobase);

            end

        end

        temp = Inf;

        for row = 1:m

            if theta(row)>0 & theta(row)

                temp = theta(row);

                outbase = XB(row);%出基变量的下标

                outrow = row;%出基变量在基变量中的位置

            end

        end

       

%%

        b(outrow) = b(outrow)/A(outrow,intobase);

        A(outrow,:) = A(outrow,:)/A(outrow,intobase);

       

        for row = 1:m

            if row ~= outrow

                b(row) = b(row) - b(outrow)*A(row,intobase);

                A(row,:) = A(row,:) - A(outrow,:)*A(row,intobase);

               

            end

        end

       

%%

        for row = 1:m

            if XB(row) == outbase;

                XB(row) = intobase;

            end   

        end

 

    end

end

%%

 

if flag == 0

    disp('LP converged to a solution x!');

    x;

    fval;

end

if flag == 1

    disp('Inf feasible solutions!');

    disp('one of the solutions is:');

end

if flag == 2

    disp('LP is unbounded!');

end

if flag == 3

    disp('No feasible point was found!');

end