二阶段法线性规划代码分享
%============================================================%');
% 二阶段法线性规划 %');
%============================================================%');
%-----------------------------------------------------------------
%说明：字母A为约束矩阵，b为资源矩阵，c（或c1）为目标函数价值系数，%
%CB为基变量系数矩阵，N为基变量下标值矩阵，t（或t1）为检验数矩阵%
%该程序假设所有决策变量都非负
%-----------------------------------------------------------------
function Twos_stage(zflag,b,A,Aflag,c)
clc
format rat
fprintf('\n');
disp('%============================================================%');
disp('% 二阶段法线性规划 %');
disp('%============================================================%');
%zflag=input('请输入所求目标函数的状态(求最大值为1，求最小值为-1):');
%b=input('请输入资源矩阵b:');
%A=input('请输入约束矩阵A:');
%Aflag=input('请输入约束方程的状态矩阵(小于号为-1，等于为0，大于号为1);');
%c=input('请输入目标函数价值系数矩阵c:');
zflag=-1;
b=[300 17.5 24];
A=[20 80;1.25 2.5;2 6];
c=[50 130];
Aflag=[-1 -1 1];
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
d=0;
continue
end
end
%-------------------------------------------------|
%以下为加入人工变量后的基变量的下标值矩阵N，约束矩
%阵A，及加入人工变量的个数n2
%-------------------------------------------------|
n2=length(find(N==0));
if n2>0 % 如果n2>0，即含有人工变量，采用二阶段法求解
for i=1:n2
A(d+3,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%daigai
CB(i)=c1(N(i)) ; %确定基变量的系数
end
[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('标准化后初始单纯形表为：');
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='次后单纯形表为：';
while k>0 & Isend==0
if (times==0)
fprintf('标准化后初始单纯形表为：');
else
fprintf('%s%d%s',temp1,times,temp2);
end
times=times+1;
[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
p=[A,b']
t,c1
%--------------------------------------------------
[m1,j]=max(t); %确定导入基，导入基的列下标为col
k=max(t);
time=0;
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
time=time+1;
end
end
if time==mA
fprintf('此问题的无最优解\n')
Isend=1
else
[m2,l]=min(f);
result1=find(f==m2);
if length(result1)~=1
[l,tmin]=min(result1)
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 %endif
end %end while
Z=c1*(X');
if Z~=0
fprintf('该问题无可行解。\n')
else
fprintf('第一阶段最优解为：\n')
find(CB==1)
~isempty(find(CB==1))
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
%******************************************************第二阶段函数完毕
%============================================================%
% 单纯形法第二阶段迭代函数 %
%============================================================%
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=(-1)*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);
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;
end
else
Isend=1;
X=zeros(1,nA);
for i=1:length(N)
X(N(i))=b(i);
end
fprintf('此问题的最优解为：\n')
c
fprintf('目标函数的最优值为：\n')
Z=c*(X')
disp('***********************************************************');
end %end if
end%end while
%*******************************************************第二阶段迭代函数结束
return