单纯形法详解及MATLAB实现，对偶单纯形法详解及MATLAB实现

单纯形法

我们以这样一个方程组做为例子，来看一下单纯形法是如何解题的
这是一个已经化成标准形式的方程组，x4和x5是我们加入的松弛变量

第一步确定一个基可行解

我们将上面方程组写作Ax=b的形式，然后从A中找出一个单位矩阵，这个单位矩阵就先做为我们的初始基矩阵，对应的解也就是基可行解。我们直接选取松弛变量的系数作为单位矩阵（这个是一定满足的），如图A矩阵对应4和5位置上的列向量组成单位矩阵
接下来很容易算出了基可行解

第二步最优性检验

也就是看一下这个解是不是最优解
简单来说，换一个x取值会不会使结果变大，如果会那这个解就不是最优解，如果不是，那我们就已经找到最优解啦。
这里当然不会是直接用眼睛看，我们是有公式的哦。
当当当，我们的检验数来啦。
首先选取基变量对应位置的c数组值
比如这里我就选取4,5位置上的c数组值，也就是0,0作为一个新的数组
我们可以记为cb=[0 0]
接下来我们对每一个变量x,xi(i=1,2,3,4,5）
计算
作为该变量所对应的检验数，m为基变量个数也就是cb数组长度，如果所有的检验数都小于等于0，那我们就找到最优解了，否则我们就要进行下一步了。

迭代运算

在开始迭代之前，我们要先看一下是不是有无界解。如果对于一个大于0的检验数，它所对应的A矩阵中的列向量元素全部小于等于0，那么我们就知道是有无界解。
否则我们就要开始迭代
具体流程我们可以用一个流程图来看一看

代码部分

首先是函数文件

function [x_opt,fx_opt,iter] = Simplex_eye(A,b,c)
%x_opt为最优解，fx_opt为最优函数值，iter为迭代次数
iter=0;%初始化迭代次数
[m,n]=size(A);%A矩阵大小
r=nchoosek(1:n,m);%选择排列
I=eye(m,m);%设置一个m阶单位矩阵，用于之后计算的比较和计算
len=length(r);
for i=1:len %从A中寻找一个单位矩阵，也就是基矩阵
    if A(:,[r(i,:)])==I
        bs=r(i,:);
        break;
    end
end
s=[1:n];
t=setdiff(s,bs);%setdiff可以计算出s数组中有，而bs数组中没有的元素
flag=0;%flag=1唯一最优解，2无穷多最优解，3有无界解
x=[];%基变量数组
while flag==0 %开始迭代
    iter=iter+1;%迭代次数加1
    x(t)=0;
    x(bs)=b;
    cb=c(bs);
    c_z=zeros(1,n); %计算检验数cj-zj
    for i=1:n
        z=sum(cb'.*A(:,i));
        c_z(i)=c(:,i)-z;
    end
    disp("----------------------第"+iter+"次---------------------------");
    All=[cb',bs',b,A];
    disp(All);
    disp(c_z);
    if all(c_z <= 0)%最优解
        x_opt=x;
        fx_opt=sum(c.*x_opt);
        if all(c_z(t) < 0)%非基变量都小于0
            flag=1;%唯一最优解
        else
            flag=2;%无穷多最优解
        end
        break;
    end
    
    [~,n1]=max(c_z);%找到最大的检验数所在位置
    disp(n1);
    disp("--------------");
    if all(A(:,n1) <= 0)%判断无界解,否则继续迭代
        x_opt=[];
        flag=3;
        break;
    end
    b1 = b./ A(:,n1);
    b1(b1<=0)=inf;%将小于等于0的数设为无穷大，
    [~,m1]=min(b1);%选出非负中的最小值对应的变量换出
    bs(m1)=n1;%n1对应为换入变量，m1对应换出变量
    t=setdiff(s,bs);
    A(:,t) = inv(A(:,bs))*A(:,t); %基矩阵的逆乘以非基矩阵
    b = inv(A(:,bs))*b;  %基矩阵的逆乘以b
    A(:,bs) = I;  %基矩阵更新为单位矩阵
end
if flag==1
    disp('唯一最优解');
    return
elseif flag==2
    disp('无穷多最优解');
    return
elseif flag==3
    disp('有无界解');
    fx_opt=inf;
    return
end

接下来是测试的脚本文件

A=[2 -3 2 1 0;1/3 1 5 0 1];
b=[15;20];
c=[1 2 3 0 0];
[x_opt,fx_opt,iter] = Simplex_eye(A,b,c)

运行结果

对偶单纯形法

函数接口文件

// DSimplex_eye.m
function [x_opt,fx_opt,iter] = DSimplex_eye(A,b,c)
% 输入参数: c为目标函数系数, A为约束方程组系数矩阵, b为约束方程组常数项
% 输出参数: x_opt最优解, fx_opt最优目标函数值, iter迭代次数
iter=0;%初始化迭代次数
[m,n]=size(A);%A矩阵大小
r=nchoosek(1:n,m);%选择排列
I=eye(m,m);%设置一个m阶单位矩阵，用于之后计算的比较和计算
len=length(r);
for i=1:len %从A中寻找一个单位矩阵，也就是基矩阵
    if A(:,[r(i,:)])==I
        ind_B=r(i,:);
        break;
    end
end
ind_N = setdiff(1:n, ind_B); 
x=[];%基变量数组
while true         % 迭代
    x(ind_N)=0;
    x(ind_B) = b;
    cB = c(ind_B);                %计算cB
    Sigma = zeros(1,n); %检验数数组
    for i=1:n
        z=sum(cB'.*A(:,i));
        Sigma(i)=c(:,i)-z;
    end
    [~,q] = min(b);               %选出最小的b，换出
    r = ind_B(q);        
    Theta = Sigma ./ A(q,:);      %计算θ
    Theta(Theta<=0) = inf;
    [~,s] = min(Theta);         %确定进基变量索引s, 主元为A(q,s)
    vals = [cB',ind_B',b,A];
    vals = [vals; NaN, NaN, NaN, Sigma];
    disp("-----------------------第"+iter+"次--------------------------")
    disp(vals);
    if all(b >= 0)         %最优解       
        x_opt = x;
        fx_opt = sum(c .* x_opt);
        return
    end
    iter=iter+1;
    % 换基
    ind_B(ind_B == r) = s;      %新的基变量索引
    ind_N = setdiff(1:n, ind_B); %非基变量索引
    % 更新A和b
    A(:,ind_N) = inv(A(:,ind_B))* A(:,ind_N);
    b = inv(A(:,ind_B))* b;
    A(:,ind_B) = I;
end

函数脚本文件

// DSimplex_eye_script.m
A=[-1 -2 -1 1 0;
    -2 1 -3 0 1];
b=[-3 -4]';
c=[-2 -3 -4 0 0];
[x_opt,fx_opt,iter] = DSimplex_eye(A,b,c)

运行结果