danchunxingfa.rar_matlab_whygpg_单纯形 matlab_单纯形法_单纯形法MATLAB

%学院：计算机学院 学号：20080353 姓名：宋勇利 %============================================================%'); % 单纯形法及其灵敏度分析程序 %'); %============================================================%'); %----------------------------------------------------------------- %说明：字母A为约束矩阵，b为资源矩阵，c（或c1）为目标函数价值系数，% %CB为基变量系数矩阵，N为基变量下标值矩阵，t（或t1）为检验数矩阵% %该程序假设所有决策变量都非负 %----------------------------------------------------------------- function danchunxingfa(zflag,b,A,Aflag,c) fprintf('\n'); disp('%============================================================%'); disp('% 单纯形法及其灵敏度分析程序 %'); disp('%============================================================%'); zflag=input('请输入所求目标函数的状态(求最大值为1，求最小值为-1):'); b=input('请输入资源矩阵b:'); A=input('请输入约束矩阵A:'); Aflag=input('请输入约束方程的状态矩阵(小于号为-1，等于为0，大于号为1);'); c=input('请输入目标函数价值系数矩阵c:'); disp('***********************************************************'); %-------------------------------------------------标准化 if zflag==-1 c=-1*c; end nb=size(b); [mA,nA]=size(A); for i=1:nb if b(i)<0 b(i)=-1*b(i); A(i,:)=-1*A(i,:); end end for i=1:mA if Aflag(i)==-1 %小于号时的处理过程 A(i,nA+1)=1; c(nA+1)=0; nA=nA+1; else if Aflag(i)==1 %大于号时的处理过程 A(i,nA+1)=-1; c(nA+1)=0; nA=nA+1; end end end %-------------------------------------------------------- %-------------------------------------------------------- %以下部分初始化基变量的下标值矩阵N,基变量的系数矩阵CB %-------------------------------------------------------- N=zeros(1,mA); CB=zeros(1,mA); v=1; for i=1:nA a=A(:,i); a1=a'; n=find(a1==0); n1=find(a1==1); a2=a1(n1); if (length(n)==mA-1)&(a2==1) %找出基变量，将下标值写入基变量的下标值矩阵N N(v)=i; v=v+1; d=n1; else continue end end %-------------------------------------------------| %以下为加入人工变量后的基变量的下标值矩阵N，约束矩 %阵A，及加入人工变量的个数n2 %-------------------------------------------------| n2=length(find(N==0)); if n2>0 % 如果n2>0，即含有人工变量，采用二阶段法求解 for i=1:n2 A(d+i,end+1)=1; %添加人工变量的约束系数 N(v)=nA+i; %将人工变量的下标写入N v=v+1; end [mA,nA]=size(A); c1=zeros(1,nA); %构造第一阶段目标函数系数矩阵 c1(nA-n2+1:nA)=-1; %人工变量的系数设为-1 for i=1:mA CB(i)=c1(N(i)) ; %确定基变量的系数 end fprintf('标准化后初始单纯形表为：'); A b c CB N %---------------------调用二阶段法单纯形法函数 [A,b,c1,N,CB]=first_stage(A,b,c1,N,CB); [A,b,c,N,CB]=second_stage(A,b,c,c1,N,CB); else %----------------如果不含人工变量，直接用单纯形法解答 for i=1:mA CB(i)=c(N(i)); end fprintf('标准化后初始单纯形表为：'); A b c CB N disp('***********************************************************'); [A,b,c,N,CB,t1]=second_stage_diedai(A,b,c,N,CB);%只调用第二阶段迭代过程 end %************************************************************* %============================================================% % 单纯形法第一阶段函数 % %============================================================% function [A,b,c1,N,CB]=first_stage(A,b,c1,N,CB) k=1; while k>0 [mA,nA]=size(A); for i=1:mA e(i)=0; i=i+1; end for j=1:nA x(j)=0; a(j)=0; t(j)=0; j=j+1; end %------------------------------------------------ for j=1:nA %start for for i=1:mA e(i)=CB(i)*A(i,j); a(j)=a(j)+e(i); i=i+1; end t(j)=c1(j)-a(j); j=j+1; end %End for %-------------------------------------------------- [m1,j]=max(t); %确定导入基，导入基的列下标为col k=max(t); if m1<=0 X=zeros(1,nA); for i=1:length(N) X(N(i))=b(i); end else %///////////////////////////////// for i=1:mA if A(i,j)>0 f(i)=b(i)/A(i,j); else f(i)=10000; %如果相除是负值，把值设为10000 end end %end for [m2,l]=min(f); w=find(CB==-1); for i=1:mA if isempty(w) break else l=min(w); end end CB(l)=c1(j);% 确定导出基 导出基的行标为 l N(l)=j; n=b(l); for i=1:mA % start for cycle ,this cycle to change the b[] if(i==l) % start if b(i)=b(l)/A(l,j); else b(i)=b(i)-A(i,j)*(n/A(l,j)); end % end the if switch i=i+1; end % end the for cycle,已改变了b 数组的值 %/////////////////////////////////////////////////////// for i=1:mA for g=1:nA P=A(i,g); if(i==l) B(i,g)=P/A(l,j); else h=A(i,j); B(i,g)=P-A(l,g)*h/A(l,j); end g=g+1; end i=i+1; end A=B; end end end fprintf('第一阶段的最终单纯形表为：'); A b c1 CB N t Z=c1*(X'); if Z~=0 fprintf('该问题无可行解。\n') else fprintf('第一阶段最优解为：\n') X end disp('***********************************************************'); return %*****************************************************第一阶段完毕 %============================================================% % 单纯形法第二阶段函数 % %============================================================% function [A,b,c,N,CB]=second_stage(A,b,c,c1,N,CB) if isempty(find(CB==-1))%判断最优解的基变量中是否含有非零的人工变量 n=length(find(c1==-1));%统计人工变量的个数 [mA,nA]=size(A); A=A(:,1:nA-n);%除去人工变量后的矩阵A [A,b,c,N,CB,t1]=second_stage_diedai(A,b,c,N,CB)%调用第二阶段迭代函数 else fprintf('此问题无解！') end return %******************************************************第二阶段函数完毕 %============================================================% % 单纯形法第二阶段迭代函数 % %============================================================% function [A,b,c,N,CB,t1]=second_stage_diedai(A,b,c,N,CB) [mA,nA]=size(A); for i=1:length(N) CB(i)=c(N(i)); end k1=1; while k1>0 for j=1:nA%(nA-n) t1(j)=c(j)-CB*A(:,j);%计算检验数 end [m3,column]=max(t1);%找出检验数中最大的值及其所在的列 k1=m3; if m3>0 B2=A(:,column); q=find(B2<0); for i=1:length(q) B2(q(i))=0; end [m4,row]=min(b./B2'); warning off MATLAB:divideByZero for i=1:mA if i==row b1(i)=b(row)/A(row,column); else b1(i)=b(i)-b(row)*A(i,column)/A(row,column); end end for i=1:mA for j=1:nA if i==row A1(i,j)=A(row,j)/A(row,column); else A1(i,j)=A(i,j)-A(row,j)*A(i,column)/A(row,column); end end end b=b1; A=A1; CB(row)=c(column); N(row)=column; else X=zeros(1,nA); for i=1:length(N) X(N(i))=b(i); end end end%end while fprintf('最终单纯形表为：'); A b c CB N t1 fprintf('此问题的最优解为：\n') X fprintf('目标函数的最优值为：\n') Z=c*(X') disp('***********************************************************'); %*******************************************************第二阶段迭代函数结束 %====================================================== % 灵敏度分析 % %======================