twolp.m

求解线性规划的两阶段法（matlab程序）里面附有相应的例子

function[z,x,p]=TwoLp(A,b,c,T)
%两阶段法求解线性规划
%入口参数A,b,c分别表示技术系数矩阵，资源限量向量，价值系数向量，入口参数T是一个向量
%向量b的维数必须等于矩阵A的行数m,向量c的维数必须等于矩阵Ａ的列数ｎ
%向量T的1至m个分量分别表示m个约束条件，0表示=，1表示<=,2表示>=
%向量T的第m+1个分量，0表示min，1表示max
%出口参数第一个记录最优值z,第二个记录最优解x
%参考例子：[z,x]=TwoLp([2,3;0,4],[26,13],[1,3],[2,2,0]);
[m0,n0]=size(A);%获得矩阵的行列数
z=0;x=[];ysn=n0;p=0;
if(m0==0)return;end
[m,n]=size(c);
if(n==1)c=c';n=m;m=1;end
if(n~=n0|m~=1)fprintf('\n错误：向量c(第3个入口参数)的维数不等于矩阵A(第1个入口参数)的列数\n\n');p=3;return;end
[m,n]=size(b);
if(m==1)b=b';m=n;n=1;end
if(m~=m0|n~=1)fprintf('\n错误：向量b(第2个入口参数)的维数不等于矩阵A(第1个入口参数)的行数\n\n');p=3;return;end
[m,n]=size(T);
if(m==1);m=n;n=1;end
if(m<=m0|n~=1)fprintf('\n错误：向量T(第4个入口参数)的维数<m+1\n\n');p=3;return;end
T=round(T);
if(m>m0)
    for(i=1:m0)if(T(i)<0|T(i)>2)fprintf('\n错误：向量T(第4个入口参数)的第%d个分量不是0，1或2\n\n',i);p=3;return;end;end
    if(T(m0+1)<0|T(m0+1)>1)fprintf('\n错误：向量T(第4个入口参数)的第%d个分量不是0，1\n\n',m0+1);p=3;return;end
end
m=m0;n=n0;wcx=0.000000001;
Control=T(m+1);
if(Control)c=-c;end
ysc=c;Vb=[];%Vb记录即基变量的序号
for(i=1:m)%将右端常数变成为非负
    if(abs(b(i))<wcx)b(i)=0;
        if(T(i)==2)A(i,:)=-A(i,:);T(i)=1;end
    elseif(b(i)<0)%不等式反向
        A(i,:)=-A(i,:);b(i)=-b(i);
        if(T(i)==2)T(i)=2;end
    end
    if(T(i))n=n+1;%化为标准型
        A(:,n)=0;
        if(T(i)==1)A(i,n)=1;Vb(i)=n;
            elseA(i,n)=-1;end
    end
end
if(n>ysn)c(n)=0;end
c0=c;c=zeros(1,n);
n1=n;%非人工变量数
for(i=1:m0)if(T(i)~=1)n=n+1;A(i,n)=1;Vb(i)=n;c(n)=1;end;end%增加人工变量A(:,n)=0;
[i,j]=size(b);
if(i==1)b=b';end
if(n==n1)c=c0;end
B=eye(m);ddbs=0;%记录迭代步数
fprintf('两阶段法求解线性规划：\n');
cB=[];for(j=1:m)cB(j)=c(Vb(j));end
while(1)ddbs=ddbs+1;%换基迭代
    bb=B*b;z=cB*bb;%计算检验数行
    s=0;cBB=cB*B;
    for(j=1:n1)%搜索正检验数
        cj=cBB*A(:,j)-c(j);
        if(abs(round(cj)-cj)<0.000001)cj=round(cj);end
        if(cj>wcx)s=j;break;end
    end
    if(s)r=0;p=B*A(:,s);minfds=Inf;%假如有正检验数
        for(i=1:m)if(p(i)>wcx&bb(i)/p(i)<minfds)minfds=bb(i)/p(i);r=i;end;end%搜索是否有轴心项
        if(r)%假如有轴心项，换基
            cB(r)=c(s);Vb(r)=s;B0=eye(m);
            for(i=1:m)B0(i,r)=-p(i)/p(r);end
            B0(r,r)=1/p(r);B=B0*B;%换基
        else p=2;fprintf('迭代步数%d，该线性规划问题没有最优解',ddbs);break;end%没有轴心项，中止计算
    else%假如没有正检验数
        if(n==n1)p=0;z=(1-2*Control)*z;fprintf('迭代步数%d,该线性规划问题的最优值为%g，最优解如下\n',ddbs,z);
                x=zeros(1,ysn);
                for(i=1:m)if(Vb(i)<=ysn)if(abs(bb(i)-round(bb(i)))<0.00001)
                            bb(i)=round(bb(i));end;x(Vb(i))=bb(i);end;end
                for(j=1:ysn)if(x(j))fprintf('    x%-4d=%g\n',j,x(j));end;end
                break%中止计算
            elseif(z>0.0001)p=1;fprintf('迭代步数%d,该线性规划问题没有可行解',ddbs);break%中止计算
            else r=0;
                for(i=1:m)if(Vb(i)>n1)r=i;break;end;end%搜索基变量中是否有人工变量
                if(r)Ar=B(r,:)*A;s=0;
                    for(j=1:n1)if(abs(Ar(j))>wcx)s=j;break;end;end
                    if(s)cB(r)=c(s);Vb(r)=s;B0=eye(m);p=B*A(:,s);
                        for(i=1:m)B0(i,r)=-p(i)/p(r);end
                        B0(r,r)=1/p(r);B=B0*B;%换基
                    else AA=[];B0=[];bb=[];
                        for(i=1:m)if(i~=r)s=s+1;AA(s,:)=A(i,:);B0(s,:)=B(i,:);bb(s)=b(i);end;end
                        clear A B T b;b=bb';
                        for(j=1:m)if(Vb(j)>Vb(r));Vb(j)=Vb(j)-1;end;end
                        s=0;for(j=1:m)if(j~=r)s=s+1;Vb0(s)=Vb(j);B(:,s)=B0(:,j);end;end
                        s=0;for(j=1:n)if(j~=Vb(r))s=s+1;A(:,s)=AA(:,j);cc(s)=c(j);end;end
                        c=cc;Vb=Vb0;clear AA B0 bb cc Vb0;
                        m=m-1;n=n-1;cB=[];for(j=i:m)cB(j)=c(Vb(j));end
                    end
                    clear Ar;
                else AA=[];
                    for(j=1:n1)AA(:,j)=A(:,j);end
                    c=c0;A=AA;clear c0 AA;n=n1;p=1;
                    cB=[];for(j=1:m)cB(j)=c(Vb(j));end
                end
            end
        end
    end