2020-11-16

最新推荐文章于 2023-05-23 14:43:37 发布

suxingfeng678

最新推荐文章于 2023-05-23 14:43:37 发布

阅读量736

点赞数

分类专栏：运筹学与最优化文章标签： matlab 线性规划

本文链接：https://blog.csdn.net/suxingfeng678/article/details/109727536

版权

运筹学与最优化专栏收录该内容

1 篇文章 0 订阅

订阅专栏

#单纯形法及大M法MATLAB实现

单纯形法是解线性规划的重要方法，matlab是强大的数学工具，非常适合用来实现单纯形法

一.单纯形法的思想

首先通过引入松弛变量与剩余变量把函数变为以下标准形式
单纯形法标准形式，其中b>=0
二.单纯形法的算法如下
（1）输入标准的系数矩阵，列向量值，目标函数系数矩阵。

（2）判断目标函数系数是否全部大于等于零，是则输出最优解与最优值，程序结束

（3）选定比值最小的主元素，进行初等行变换，直到出现单位矩阵与基变量的检验数为零则结束此步骤。

回到（2）

三.大M法

大M法是为解决b非负与单位子块不可兼得的情况。即把系数矩阵化为标准形式不管b列正负。大M法算法步骤与单纯形法一样，在单纯形法加入M乘以单位矩阵，只需要输入时换为大M法的格式，其余步骤与第一题单纯形法一模一样。设M为比其他值都大一千倍以上的正数，目标函数系数也写入M。

四.单纯形法的特殊情况
（1）当使用大M法求极大值目标函数时，b列元素有为零的情况，则无解。
如果是单纯形法，因为引入了一个单位子块，初始状态时一定是有解的。
（2）若基本可行解是非退化的，且c列中所有检验数大于等于零并且至少有一个等于0，那么有无穷个最优解。
（3）单纯形表中有一列全为<=0的，线性规划无界。

matlab代码实现如下

function[A,m,flag]=simplex(A,m,b,c)%A为标准系数矩阵
A=[3 5 0 1 0 0 0;0 1 0 0 1 0 0;8 5 0 0 0 1 0;1 1 -1 0 0 0 1];
b=[150;20;300;50];
M=100000;        %大M法
c=[-50 -40 0 M M M M]
[a1 a2]=size(A)    %取A的行数a1,列数a2；
m=[((a2-a1)+1):a2] %列出基变量的下标号；
n=[1:(a2-a1)]         %列出非基变量下标；
A1=[A b]
flag=1;
k=0;
rl=0;
r=1;
while flag
b=A1(:,a2+1)
for i=1:a1
    if m(i)==rl
        m(i)=k; 
    end
end
for i=1:(a2-a1)
    if n(i)==k
        n(i)=rl;
    end
end
E=A1(:,m)
%建立初始单纯形表；
cb=c(m(1));
xb=(inv(E))*b
for i=2:a1
    cb=[cb c(m(i))]%取基向量cb;
end
z0=cb*xb          %计算目标函数,即基本可行解；
rr=zeros(1,a2)
%计算检验数；
for j=1:(a2-a1)                          %计算检验数；
    cy=0;
    for i=1:a1
        cy=cy+c(m(i))*A1(i,n(j));
    end
    z(n(j))=cy
    rr(n(j))=c(n(j))-z(n(j))
end
A2=[rr -z0]
AA=[A1;A2]    %生成初始单纯形表；
n
if (min(rr))>=0    %判断是否有解。是否是一个最优解，如果是则显示变量取值x以及最优值z；
    for q=1:(a2-a1)
         if(rr(n(q))==0&&~all(A1(:,a2)))
             disp('有无穷多个最优解,下面为其中之一')
        end  

   end
    z=z0
    x(m)=xb
    x=x(:)
    disp('变量最优解为：x*=')
    disp(x)
    disp('最优解值为：z*=')
    disp(z)
    flag=0;
   break;
end
%决定进基矢量ak;

%决定离基矢量ar和主元素yrk;
k=min(find(rr==min(rr(find(rr<0)))))
ak=A1(:,k)
if max(A1(:,k))<=0        %如果有一列全部 小于等于0的值则目标函数是无界的，无最优解；
    disp('没有边界')
    flag=0;
    break;
end
%if max(A2(:,k))<=0        %如果有一列（包括c列)全部小于等于0的值则目标函数是无可行解；
%    disp('无可行解')
%    flag=0;
%    break;
%end
for i=1:a1
    NN(i)=(A1(i,a2+1))/(A1(i,k))
end
r=find(NN==min(NN(find(NN>0))))     %找主元素，即NN中最小的那个对应的列
    yrk=A1(r,k)            %主元素为yrk； 
    rl=m(r)                %离基矢量下标为rl；
   % ar=A1(:,rl)          %离基矢量为ar；
%行变换，参考了数值分析的高斯消去法
for i=1:a1
    if i~=r
        A1(i,:)=A1(i,:)-A1(r,:)*A1(i,k)/A1(r,k)
    else
        A1(r,:)=A1(r,:)/yrk;
    end
end
A1
end