MMSimplex.rar_MMSimplex_MMSimplex.m文件_单纯形法_单纯形法MATLAB_运筹学

function [xstar,fxstar,iter,A0,IB,IN,SA,xSA,InverseofB,exitflag]=MMSimplex(A,b,c) A0=A; %系数矩阵 ATTENTION : 该矩阵已经经过了初步的调整——满足标准化形式 ；采用两阶段法 [m,n]=size(A0); E=eye(m); %需要解释 IB=zeros(1,m); %最优基变量下标 SA1=zeros(1,n); %松弛变量或者剩余变量下标 IR1=zeros(1,m); %需要在本code的后面用到 IR=1:m; k=0; tic; %star stopwatch timer to measure performance for i=1:m for j=1:n if A0(:,j)==E(:,i) %判断系数矩阵的第 j 列是否为一单位向量 IB(i)=j; %确定基变量下标 IR1(i)=i; %记录原来系数矩阵中已有的E（：，i） SA1(i)=j; %松弛变量 elseif A(:,j)==(-E(:,i)) SA1(i)=j; %剩余变量 end end end s1=find(SA1~=0); %注意不等于号；find函数(举例子)：用于返回元素的所在 位置 if length(s1)~=0 %有剩余变量有松弛变量 for i=1:length(s1) SA(i)=SA1(s1(i)); %记录松弛变量和剩余变量的下标 ; 存在的意义SA与SA1的列数不同___去“杂质” end else SA=[]; end IR=find(IR~=IR1); %IR=1：m； 系数矩阵中 不是 单位向量的列——需要通过人工变量补齐的length（s）个人工变量 s=find(IR~=0); %length(s)是否=length(IR)? for p=1:length(s) A0(:,n+p)=E(:,IR(p)); %添加人工变量——缺少哪个补充哪个 IB(IR(p))=n+p; %更换下标 IR(p)=n+p; %需要添加的第p个人工变量的下标记录为p+n——从此以后，IR 记录人工变量的下标 end A0=[b,A0]; %合并b列和加入了人工变量之后的系数矩阵 flag=0; %通俗的理解成一张“通行证”，理解成一种标志 while(length(IR)~=0)&(flag==0) c0=zeros(n+length(s),1); %两阶段法，并且要把第一阶段的目标函数也化成标准型，从而和后面的程序对应（检验数） c0(IR)=-ones(length(s),1); %此处IR表示人工变量的下标; 该语句表示生成仅有人工变量的系数为-1的目标函数（注意min z与max z的区别） N=1:n+length(s); %表示所有变量（包括原有变量、松弛变量、剩余变量、人工变量）的下标 N(IB)=[]; %IB现阶段已经表示所有基变量的下标；空数组可以用于数组的删除，删去所有基变量的下标。找出所有非基变量的下标 IN=N; %删除了基变量下标的数组就是非基变量下标组成的矩阵 IN(find(IN==0))=[]; x(IN)=zeros(1,length(IN)); %令初始基可行解里的所有非基变量为0，所有基变量等于对应的b值 x(IB)=A0(:,1)'; cB=c0(IB); %记录基变量对应的价值系数矩阵为cB sigma=c0'-cB'*A0(:,2:n+length(s)+1); %A0第一列为b； sigma是一个行向量 t=length(find(sigma>0)); while t~=0 [sigmaJ,jj]=max(sigma); %表示位置——列；确定在sigma这一行向量中max sigma的位置（即A0中的jj+1列） tt=find(A0(:,jj+1)>0); %找到该变量对应的系数大于0的行 kk=length(tt); if kk==0 %课本P32 disp('原问题无解1'); xstar=[];fxstar=[];A0=[];IB=[];iter=k; flag=1; else %证明存在非基变量符合换入变量的条件，下面开始寻找符合条件的换出变量 theta=zeros(1,kk); %为需要计算的theta开辟向量存储空间 for i=1:kk theta(i)=A0(tt(i),1)/A0(tt(i),jj+1); end [thetaI,ii]=min(theta); %确定在theta这一行向量中min theta的位置 Temp=tt(ii); %表示位置——行；确定换出变量 for i=1:m if i~=Temp A0(i,:)=A0(i,:)-(A0(Temp,:)/A0(Temp,jj+1))*A0(i,jj+1); %《MATLAB应用P40》 else A0(Temp,:)=A0(Temp,:)/A0(Temp,jj+1); end end %以上为旋转运算 TT=IB(Temp); %TT为中间变量；数组IB中第Temp个IB 《MATLAB应用P43》； TT表示被换出量的下标 IB(Temp)=jj; %非基变量换入 for i=1:length(IR) if IR(i)==TT %如果已经把第i个人工变量从基中换出的话 IR(i)=0; %那么该人工变量不再有影响，该位置赋值为0 end end d=find(IR==0);IR(d)=[]; %这里记录的是人工变量的变化：上面的for语句和该语句共同表达：换出人工变量并删除 IN(jj)=TT; %这里是迭代中形成新的非基变量组合的过程；原来的基变量（被换出的）下标成为了现在的非基变量下标 x(IB)=A0(:,1)'; IN(find(IN==0))=[]; x(IN)=zeros(1,length(IN));cB=c0(IB); %新的基可行解及价值系数 sigma=c0'-cB'*A0(:,2:n+length(s)+1); t=length(find(sigma>0)); %再次计算检验数并假设检验中有t个大于零的检验数 end k=k+1; %即已经完成第一次迭代 end %完成第51行开始的while循环，所有的迭代步骤完成；即不再满足第51行的循环条件：sigma<0 if sum(x(IR))~=0 %两阶段法第一阶段的判断指标 disp('原问题无解2'); %此时没有检验数小于零，但第一阶段有最优解，从而原问题无解 xstar=[];fxstar=[];A0=[];IB=[];iter=k; flag=2;exitflag=flag; else x=x(1:n); end end if (flag==1)|(flag==2) %flag=1出现在58行，flag=2出现在95行 return else IB; N=1:n; N(IB)=[ ]; IN=N; IN(find(IN==0))=[ ]; x(IN)=zeros(1,length(IN)); cB=c(IB); A0=A0(:,1:n+1); sigma=c'-cB'*A0(:,2:n+1); t=length(find(sigma>0)); while (t~=0)&&(flag==0) [sigmaJ,jj]=max(sigma); tt=find(A0(:,jj+1)>0); kk=length(tt); if kk==0 disp('原问题为无界解3'); xstar=[]; fxstar=[]; A0=[]; IB=[]; iter=k; flag=1; else theta =zeros(1,kk); for i=1:kk theta(i)=A0(tt(i),1)/A0(tt(i),jj+1); end [thetaI,ii]=min(theta);Temp=tt(ii); for i=1:m if i~=Temp A0(i,:)=A0(i,:)-(A0(Temp,:)/A0(Temp,jj+1))*A0(i,jj+1); else A0(Temp,:)=A0(Temp,:)/A0(Temp,jj+1); end end TT=IB(Temp); IB(Temp)=jj; IN(jj)=TT; x(IB)=A0(:,1)'; N=1:n; N(IB)=[]; IN=N; IN(find(IN==0))=[]; x(IN)=zeros(1,length(IN)); cB=c(IB); sigma=c'-cB'*A0(:,2:n+1); t=length(find(sigma>0)); end k=k+1; end end if flag==1 xstar=[]; fxstar=[]; A0=[]; IB=[]; iter=k; disp('原问题为无界解4');exitflag=flag; elseif flag==2 xstar=[]; fxstar=[]; A0=[]; IB=[]; iter=k; disp('原问题无解5'); exitflag=flag; else xstar=zeros(1,n); xstar(IB)=A0(:,1)'; %xstar为最优解 fxstar=xstar(IB)*c(IB); %fxstar 为最优值——求max z的过程 iter=k; B=A(:,IB); InverseofB=inv(B); xSA=x(SA); exitflag=flag; end toc;