function Two_stage(zflag,b,A,Aflag,c)
clear;clc
format rat
fprintf('\n');
%zflag=input('请输入所求目标函数的状态(求最大值为1，求最小值为-1):');
%b=input('请输入资源矩阵b:');
%A=input('请输入约束矩阵A:');
%Aflag=input('请输入约束方程的状态矩阵(小于号为-1，等于为0，大于号为1);');
%c=input('请输入目标函数价值系数矩阵c:');
zflag=1;
b=[9 15 5];
A=[5 3 1;-5 6 15;2 1 1];
c=[10 15 12];
Aflag=[-1 -1 1];
%-------------------------------------------------标准化
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
num1=find(Aflag==-1);
num2=find(Aflag==0);
num3=find(Aflag==1);
N=[];
c1=zeros(1,length(c)); %构造第一阶段目标函数系数矩阵
for i=1:length(num1)
A(num1(i),nA+1)=1;
c(nA+1)=0;
c1(nA+1)=0;
nA=nA+1;
N(i)=nA;
end
for i=1:length(num3)
A(num3(i),nA+1)=-1;
A(num3(i),nA+2)=1;
c(nA+1)=0;
c1(nA+1)=0;
c1(nA+2)=-1;
nA=nA+2;
N(end+1)=nA;
end
for i=1:length(num2)
A(num2(i),nA+1)=1;
c1(nA+1)=-1;
nA=nA+1;
N(end+1)=nA;
end
%以下部分初始化基变量的下标值矩阵N,基变量的系数矩阵CB
%--------------------------------------------------------
CB=zeros(1,mA);
[mA,nA]=size(A);
for i=1:mA
CB(i)=c1(N(i)) ; %确定基变量的系数
end
if (length(num2)+length(num3))~=0
[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
fprintf('标准化后初始单纯形表为：');
c
p=[CB',N',b',A]
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; % k是
times=0;
Isend=0;
temp1=' 第';temp2='次后单纯形表为：';
b=b';
while Isend==0
if (times==0)
fprintf('标准化后初始单纯形表为：');
else
fprintf('%s%d%s',temp1,times,temp2);
end
times=times+1;
[mA,nA]=size(A);
t=zeros(1,mA);
for j=1:nA %start for
t(j)=c1(j)-CB*A(:,j);
end %End for
c1
p=[CB',N',b,A]
t
[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
break; %跳出while 循环
else
times2=0;
for i=1:mA
if A(i,j)>0
f(i)=b(i)/A(i,j);
else
f(i)=10000; %如果相除是负值，把值设为10000
times2=times2+1;
end
end
if times2==mA
fprintf('此问题的无最优解\n')
Isend=1;
break; %跳出if m1<0
else
[m2,l]=min(f);
CB(l)=c1(j);% 确定导出基 导出基的行标为 l
N(l)=j;
B=[];
for i=1:length(N)
B=[B,A(:,N(i))];
end
b=B\b ; % end the for cycle,已改变了b 数组的值
A=B\A;
end %endif
end %endif
end %end while
Z=c1*(X');
if Z~=0
fprintf('人工变量仍在基变量中，该问题无可行解。\n')
else
fprintf('第一阶段最优解为：\n')
X
end
disp('***********************************************************');
%*****************************************************第一阶段完毕
% 单纯形法第二阶段函数 %
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
% 单纯形法第二阶段迭代函数 %
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;
times=0;
Isend=0;
while Isend==0
if (times==0)
fprintf('初始单纯形表为：');
else
fprintf('%s%d%s\n','第',times,'次迭代的单纯性表');
end
times=times+1;
for j=1:nA%(nA-n)
t1(j)=c(j)-CB*A(:,j); %计算检验数
end
p=[CB',N',b,A]
c,t1
[m3,column]=max(t1);%找出检验数中最大的值及其所在的列
k1=m3;
j=column;
ftimes=0;
if m3>0
for i=1:mA
if A(i,j)*b(i)>0
f(i)=b(i)/A(i,j);
else
f(i)=10000; %如果相除是负值，把值设为10000
ftimes=ftimes+1;
end
end
if ftimes==mA
fprintf('此问题的无最优解\n')
Isend=1;
else
[m4,row]=min(f);
CB(row)=c(column);
N(row)=column;
B=[];
for i=1:length(N)
B=[B,A(:,N(i))];
end
b=B\b; % end the for cycle,已改变了b 数组的值
A=B\A;
end
else
Isend=1;
X=zeros(1,nA);
for i=1:length(N)
X(N(i))=b(i);
fprintf('此问题的最优解为：\n')
X
fprintf('目标函数的最优值为：\n')
Z=c*(X')
disp('***********************************************************');
end
end %end if
end%end while
%*******************************************************第二阶段迭代函数结束
