不了解单纯形法的两阶段法的matlab实现？看这篇就够了

---

前言

最优化算法两阶段——单纯形法的matlab实现：
要读懂本文代码需要了解：nchoosek，setdiff，eye等函数在matlab中的作用，以及/符号在矩阵运算中的作用。

---

一、单纯形法表格

在高等教育出版社《最优化方法》一书中提到的单纯形法表格如下图所示：

其中：c为目标函数系数, A为约束方程组系数矩阵, b为约束方程组常数项。

1.1可立即读出最优解和最优值的表格具备的特点

① 中心部位有单位子块
② 右列元素非负
③ 底行相应于单位子块的位置为0
④ 底行其他元素非负

---

二、两阶段——单纯形法的步骤（流程图）

---

三、两阶段——单纯形法的matlab实现

3.1 两阶段——单纯形法matlab代码

文件1——dcxf2jd.m

function x=dcxf2jd(A,b,c)
%% 介绍
% 两阶段单纯形法求线性规划最优解，针对于A中未必刚好有一个m阶单位矩阵
% min cx   s.t Ax=b,x>=0,注意A是m*n矩阵，b是m维列向量， c是n维行向量
%% 第一阶段得到辅助问题的单纯形表
[m,n]=size(A);
basis=((n+1):(n+m))';                                                      % 获取基变量的下标
cB=ones(1,m);                                                              % 获取基变量对应的cB
zj_cj=cB*A;                                                                % 获取非单位矩阵对应检验数
f=cB*b;                                                                    % 计算当前的函数值
form=[basis A eye(m) b;0 zj_cj zeros(1,m) f];                              % 构建初始化单纯形表,form是(m+1)*(m+n+2)的矩阵,form(m+1,2:(m+n+1)) 是检验数向量
%% 进入循环while求解目标函数（所有人工变量求和）的最优解
 while any(form(m+1,2:(m+n+1))>0)                                          % m+1行，2到m+n+1列中的检验数>0,如果检验数中有>0的，表明把这个>0的分量作为进基变量，函数值还可以更小，没有找到最优值
    zj_cj=form(m+1,2:(m+n+1));                                             % 把所有的判别数的值赋给zj_cj
    [~ , k]=max(zj_cj);                                                    % 求得所有判别式中最大的那个数的下标是k,也就是进基变量的下标,从左往右进行比较，如果有两个一样的最大值，则选择最左边的那个
       yk=form(1:m,k+1);                                                   % 由于form矩阵还有第一列是基变量下标basis，所以如果在form矩阵中引用这一列，标号应该是k+1,yk位1到m行，k+1列的数
    [yk1 ,J]=find(yk'>0);                                                  % 找到yk中>0的那些数
    if isempty(yk1)                                                        % 如果没有大于0的，表明不存在最优解
        error('问题不存在最优解');                                         
    end
    b_=form(1:m,m+n+2);                                                    % 对应的是b的一个列向量
    [~, J1]=min(b_(J)./yk(J));                                             % J1表示b_(J)./yk(J)中最小数所在列位置(./表示的是对应元素的除法运算)
    r=J(J1);                                                               % 找到进基列的那个变化后非零分量的行数r
    y=form(:,k+1);                                                         % y是进基变量对应的整一列向量
    E=qiuni(y,r);                                                          % 下面要对整个单纯性表（不包括basis)做初等行变换，使得y向量变成(0,0,...1,0,...0)T,1在r个分量位置
                                                                           % 由于初等行变换相当于左乘一个可逆矩阵，这一个子函数qiuni就是求出这一个可逆矩阵
    whole1=E*form(:,2:(m+n+2));                                            % 做初等行变换（左乘可逆矩阵）
    basis(r)=k;                                                            % 更新basis基变量下标
    form=[[basis;0] whole1];                                               % 更新form                 
end
%% 第一阶段求解最优解结束，进入判断三种情况
% ①最后的最优值大于0，无可行解
% ②最后的最优解等于0，基变量全在x1,x2,...,xn中，则x=(x1,x2,...,xn)T就是初始基可行解
% ③最后的最优解等于0，基变量不全在x1,x2,...,xn中，仍然有人工变量yk是基变量，则要把yk变为非基变量
if abs(form(end,end))==0                                                   % 如果最优值=0，表明原问题存在可行解，否则不存在
[~,s2]=find(basis'>n);                                                     % s2保存basis中大于n的下标，由于大于n的基变量为人工变量，要用原有的变量将其驱赶
if ~isempty(s2)                                                            % 如果s2是空，则不存在大于n的下标，直接把form的最后一行和人工变量所在列去掉就行，见下面
    p=size(s2,2);                                                          % p是s2中的列数，在这相当于求s2元素个数，表示基变量的选取选择人工变量的有s2个
    del=zeros(1,p);                                                        % 如果有一个人工变量是基变量但是这一行对应的所有原来的变量的值均为0，即无法驱赶，则将这一行整个删掉,用del记录
    for i=1:p
        ind=find(form(s2(i),2:n+1));                                       % 寻找这一行中x1到xn列中非零元素
        if ~isempty(ind)                                                   % 如果没找到，则删除这一行，见下面，找到了，就取第一个
             k1=ind(1);                                                    % 在form中位置是k1+1;
             E1=qiuni(form(:,k1+1),s2(i));                                 % 驱赶的过程和上面进基变量替换离基变量过程类似，初等行变换，
             whole1=E1*E*form(:,2:(m+n+2));                                % 相当于左乘一个矩阵
             basis(s2(i))=k;                                               % 更新基变量下标
             form=[[basis;0] whole1];                                      % 更新form单纯形表
        else
             del(i)=1;                                                     % 记录要删除的这一行的位置
        end
    end
    ind2=setdiff(1:m,s2(logical(del)));                                    % setdiff(a,b) 相当于集合减法，a-b=a-(a交b) a=1:m而form有m+1列，相当于同时把form的最后一行去掉,ind2得到删除后剩下的下标
    form=form(ind2,[1:(n+1) m+n+2]);                                       % 新的form矩阵
else
    form=form(1:(end-1),[1:(n+1) m+n+2]);                                  % 如果不需要驱赶，则直接把form的最后一行和人工变量所在列去掉
end
else
    error('原问题不存在可行解');
end
%% 通过第一阶段找到初始基本可行解后把form矩阵的判别式列补上
[m1 , n1]=size(form);
zj_cj2=c(form(:,1))*form(:,2:(n1-1))-c;                                    % 计算判别式，注意这里把基变量和非基变量放在一起算，如果是基变量，zj-cj2这样算的结果就等于0
b_2=c(form(:,1))*form(:,n1);                                               % 计算函数值，这里选出基变量所对应的c的值c(form(:,1))来求判别式
form2=[form;[0 zj_cj2 b_2]];                                               % form2是第二阶段的单纯形表,form2 的大小: m1+1 * n1
basis=form(1:m1,1);                  
%% 第二阶段求最优解，和第一阶段类似
while any(form2(m1+1,2:(n1-1))>0)                                          % m1+1行，2到n1-1列中的检验数>0,如果检验数中有>0的，表明把这个>0的分量作为进基变量，函数值还可以更小，没有找到最优值
     zj_cj=form2(m1+1,2:(n1-1));                                           % 把所有的判别数的值赋给zj_cj
     [~ , k]=max(zj_cj);                                                   % 求得进基变量的下标,从左往右进行比较，如果有两个一样的最大值，则选择最左边的那个
    yk=form2(1:m1,k+1);                                                    % 在form2里面进基变量的位置是k+1;
    [yk1 , J]=find(yk'>0);                                                 % 注意 find函数中 行向量和列向量的返回值不同，如果是列向量，则第一个返回值是下标，行向量是第二个返回值是下标
    if isempty(yk1)                                                        % 找到yk中>0的那些数
        error('问题不存在最优解');                                          % 如果没有大于0的，表明不存在最优解
    end
    b_=form2(1:m1,n1);                                                     % 对应的是b的一个列向量
    [~, J1]=min(b_(J)./yk(J));                                             % J1表示b_(J)./yk(J)中最小数所在列位置
    r=J(J1);
    y=form2(:,k+1);                                                        % y是进基变量对应的整一列向量
    E=qiuni(y,r);                                                          % 求出可逆矩阵E
    whole3=E*form2(:,2:(n1));                                              % 做初等行变换（左乘可逆矩阵）
    basis(r)=k;                                                            % 更新basis
    form2=[[basis;0] whole3];                                              % 更新表form2
end
%% 求最优解所对应的x
[m3,n3]=size(form2);
basis=form2(1:m3-1,1);
 s=n3-2;
 x=zeros(1,s);
 for i=1:m3-1
 x(basis(i))=form2(i,n3);
 end
end

文件2——qiuni.m

function E=qiuni(y,num)
% dcxf2jd的子函数
%给出一个向量y和一个分量的位置num，求得把这个向量变成(0,0,0,...,1,0,...0)，1在num位置向量的左乘矩阵E
% y是列向量！！！
m=size(y,1);                                                               % 获取y的行数
if abs(y(num))==0                                                          % 此时y(num)=0，不存在这样的E
    E=-1;
else
    E=eye(m);
    y1=-y;
    y1(num)=1;
    y1=y1./(y(num));
    E(:,num)=y1;
end
end

---

3.2测试例题

clc,clear;
A=[1,1,1,1,0;
   -2,1,-1,0,-1;
   0,3,1,0,0];
b=[4,1,9]';
c=[3,0,-1,0,0];
x=dcxf2jd(A,b,c)
a=c.*x;
fm=sum(a(:))

3.3结果

以上就是两阶段——单纯形法的matlab实现，有什么疑惑的地方可以私信或者评论区提问，需要流程图的也可以私信我。
本文参考的相关文章链接：最优化算法单纯形法的matlab实现（单纯形法看这一篇就够了)