自编MATLAB版单纯性算法 可以列出单纯形表以及其他相关数据

function [dcxb,x,fval,exitflag,flag]=simplex(f,A,b,Aeq,beq)

%本程序相关说明：

%(1)本函数基本使用方法与MATLAB自身函数使用方法类似

%(2)dcxb表示整个过程的单纯形表，flag=1表示没有出现退化情况，flag=-1表示出现退化情况，改为依照勃兰德法则确定进基变量和出基变量

%(3)程序运行结果dcxb中inf无意义，表的第一行两inf中间的数字表示变量的下标，本程序默认变量下标从1开始

%(4)dcxb其余各行如果仅仅行首和行尾出现inf，则两inf中间内容为检验数

%(5)dcxb第一列列相邻两inf间的内容为基变量的下标

%(6)exitflag为负表示无解，此时x,fval并不是相应的基可行解和最有值，仅仅代表程序中止时的基解和该基解所对应的函数值

%(7)函数名暂时命名为simplex，simplex是谷歌翻译给出的单纯形法的英文翻译，不知是否为专业术语

%(8)flag=-1时dcxb中会出现一行全为inf，无实际意义，仅表示该行一下依照勃兰德法则确定进基变量和出基变量(参见用法示例(2))

%(9)其他参数与linprog中相应参数意义相似

%(10)exitflag=1表示有唯一最优解，exitflag=2表示有无穷解，exitflag=-1表示目标函数无下界，exitflag=-2表示可行域为空

%(11)flag=2表示用到了两阶段法，此时程序运行结果将出现“第一阶段单纯形表为”“第二阶段单纯形表为”字样，并显示两个阶段各自的单纯形表，详见用法示例(3)

%用法示例(1)：题目为束金龙编的《线性规划理论与模型应用》P44第23题(1)小题

%在命令窗口依次输入以下语句：

%f=[1 -2 1 -3]';

%A=[1 1 3 1;0 -2 1 1;0 -1 6 -1];

%b=[6 3 4]';

%[dcxb,x,fval,exitflag,flag]=simplex(f,A,b)

%用法示例(2):题目为束金龙编的《线性规划理论与模型应用》P23例1.8

%命令窗口输入的命令为：

%f=[-0.75 20 -0.5 6 0 0 0]';

%Aeq=[0.25 -8 -1 9 1 0 0;0.5 -12 -0.5 3 0 1 0;0 0 1 0 0 0 1];

%A=[];

%b=[];

%beq=[0 0 1]';

%[dcxb,x,fval,exitflag,flag]=simplex(f,A,b,Aeq,beq)

%用法示例(3):题目为课本《线性规划理论与模型应用》P29例1.9

%命令窗口输入的命令为：

%f=[1 -2]';

%A=[-1 -1;1 -1;0 1];

%b=[-2 -1 3]';

%[dcxb,x,fval,exitflag,flag]=simplex(f,A,b)

flag=1;

if nargin < 5, beq = [];

if nargin < 4, Aeq = [];

end

end

[A_m,A_n]=size(A);

[Aeq_m,Aeq_n]=size(Aeq);

bb=[b;beq];

m=A_m+Aeq_m;

n = max([length(f),A_n,Aeq_n]); % In case A is empty

if isempty(f), f=zeros(n,1); end%4

if isempty(A), A=zeros(0,n); end

if isempty(b), b=zeros(0,1); end

if isempty(Aeq), Aeq=zeros(0,n); end

if isempty(beq), beq=zeros(0,1); end

if(all(b>=0)==0)

flag=3;

b_len=length(b);

k=1;

for r=1:inf

if(k<=b_len)

if b(k)<0

Aeq=[Aeq,zeros(Aeq_m,1);A(k,:),1];%1%5

beq=[beq;b(k)];%2

A=[A(1:k-1,:);A(k+1:A_m,:)];%3

A=[A,zeros(A_m-1,1)];

b=[b(1:k-1);b(k+1:A_m)];

k=k-1;

b_len=b_len-1;

A_m=A_m-1;

Aeq_m=Aeq_m+1;

Aeq_n=Aeq_n+1;

f=[f;0];

else

Aeq=[Aeq,zeros(Aeq_m,1);A(k,:),1];

beq=[beq;b(k)];

A=[A(1:k-1,:);A(k+1:A_m,:)];

A=[A,zeros(A_m-1,1)];

b=[b(1:k-1);b(k+1:A_m)];

k=k-1;

b_len=b_len-1;

A_m=A_m-1;

Aeq_m=Aeq_m+1;

Aeq_n=Aeq_n+1;

f=[f;0];

end

k=k+1;

else break;

end

end

A=[];

b=[];

end

for k=1:Aeq_m

if beq(k)<0

Aeq(k,:)=-Aeq(k,:);

beq(k)=-beq(k);

end

end

[A_m,A_n]=size(A);

[Aeq_m,Aeq_n]=size(Aeq);

bb=[b;beq];

m=A_m+Aeq_m;

n = max([length(f),A_n,Aeq_n]);

if isempty(f), f=zeros(n,1); end

if isempty(A), A=zeros(0,n); end

if isempty(b), b=zeros(0,1); end

if isempty(Aeq), Aeq=zeros(0,n); end

if isempty(beq), beq=zeros(0,1); end

AA=[A;Aeq];

pq=0;

PQ=zeros(n+m,1);

QP=zeros(n,1);

for k=1:n

D=AA(:,k);

[maxk,maxK]=max(D);

[mink,minK]=min(D);

if(maxk~=0)%这种情况是不正确的

if(mink==0)

D(maxK)=0;

maxk=max(D);

if(maxk==0)

pq=pq+1;

PQ(maxK)=k;

QP(pq)=maxK;

for s=1:pq-1

if(A(:,s)==A(:,pq))

pq=pq-1;

end

end

end

end

end

end

PQ=PQ(1:m);

QP=QP(1:pq);

W=eye(m);

for k=1:pq

W(:,QP(k))=zeros(m,1);

end

WW=zeros(m,0);

kk=1;

for k=1:m

if(W(:,k)==zeros(m,1))

else

for r=1:m

if PQ(r)==0

PQ(r)=n+kk;

break;

end

end

WW=[WW,W(:,k)];

kk=kk+1;

end

end

AA=[AA,WW];

Ab=[AA,bb];

FF=zeros(1,m+n-pq);

for r=n+1:m+n-pq

for s=1:m

if AA(s,r)~=0

FF=FF+AA(s,:)/AA(s,r);

end

end

end

TT=f;

FF=[FF(1:n),zeros(1,m-pq)];

if flag==3

ff=-FF';

T=[f;zeros(m-pq,1)];

else

ff=[f;zeros(m-pq,1)];

end

F=ff;

B=PQ;

BB=B;

dcxb=[inf,1:m+n-pq,inf;B,Ab;inf,ff',inf];

for loop=1:inf

[minff,i]=min(ff);

if(flag~=1)

for k=1:i

if(ff(k)<0)

i=k;

break;

end

end

end

C=inf*ones(m,1);

if(minff>=0)

exitflag=1;

break;

else

if(AA(:,i)+abs(AA(:,i))==0)%注意等号

exitflag=-1;%exitflag为负表示无解

dcxb=[dcxb;inf(1,m+n-pq)];

break;

else

for h=1:m

if(AA(h,i)>0)

C(h)=bb(h)/AA(h,i);

end

end

[minC,j]=min(C);

B(j)=i;

Ab(j,:)=Ab(j,:)/Ab(j,i);

ff=ff-AA(j,:)'*ff(i)/AA(j,i);

if flag>1

T=T-AA(j,:)'*T(i)/AA(j,i);

end

for k=1:m

if(k~=j)%不等号<>不对吗？

Ab(k,:)=Ab(k,:)-Ab(j,:)*Ab(k,i)/Ab(j,i);

end

end

bb=Ab(:,n+m-pq+1);

AA=Ab(:,1:n+m-pq);

dcxb=[dcxb;B,Ab;inf,ff',inf];

BB=[BB,B];

end

end

if flag==1

for r=1:(loop)

if(BB(:,r)==BB(:,loop+1))

flag=-1;

dcxb=[dcxb;inf*ones(1,m+n+2-pq)];

break;

end

end

end

end

if flag==3

x=zeros(length(ff),1);

for k=1:m

x(B(k))=bb(k);

end

if x'*ff~=0

exitflag=-2

end

fprintf('本题采用两阶段法求解\n第一阶段单纯形表为：\n')

dcxb

flag=2;

BB=[];

ff=T(1:n);

end

if flag==2

F=ff;

AA=AA(:,1:n);

Ab=[Ab(:,1:n),bb];

BB=B;%%%%%

dcxb=[inf,1:n,inf;B,Ab;inf,ff',inf];

for loop=1:inf

[minff,i]=min(ff);

if(flag==-1)

for k=1:i

if(ff(k)<0)

i=k;

break;

end

end

end

C=inf*ones(m,1);

if(minff>=0)

exitflag=1;

break;

else

if(AA(:,i)+abs(AA(:,i))==0)%注意等号

exitflag=-1;%exitflag为负表示无解

dcxb=[dcxb;inf(1,m+n-pq)];

break;

else

for h=1:m

if(AA(h,i)>0)

C(h)=bb(h)/AA(h,i);

end

end

[minC,j]=min(C);

B(j)=i;

Ab(j,:)=Ab(j,:)/Ab(j,i);

ff=ff-AA(j,:)'*ff(i)/AA(j,i);

for k=1:m

if(k~=j)%不等号<>不对吗？

Ab(k,:)=Ab(k,:)-Ab(j,:)*Ab(k,i)/Ab(j,i);

end

end

bb=Ab(:,n+1);

AA=Ab(:,1:n);

dcxb=[dcxb;B,Ab;inf,ff',inf];

BB=[BB,B];

end

end

if flag==1

for r=1:(loop)

if(BB(:,r)==BB(:,loop+1))

flag=-1;

dcxb=[dcxb;inf*ones(1,n+2)];

break;

end

end

end

end

end

if flag==2

x=zeros(n,1);

F=TT;

fprintf('第二阶段单纯形表为：\n')

else

x=zeros(m+n-pq,1);

end

for k=1:m

x(B(k))=bb(k);

end

if sum(ff==0)~=m

exitflag=2;

end

fval=x'*F;

程序可能会有瑕疵之处  因为我并没有怎么学过程序  MATLAB也是刚接触  若有错误请见谅  但求抛砖引玉