单纯形法详解及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,:)])==Ibs=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;%迭代次数加1x(t)=0;x(bs)=b;cb=c(bs);c_z=zeros(1,n); %计算检验数cj-zjfor i=1:nz=sum(cb'.*A(:,i));c_z(i)=c(:,i)-z;enddisp("----------------------第"+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)%非基变量都小于0flag=1;%唯一最优解elseflag=2;%无穷多最优解endbreak;end[~,n1]=max(c_z);%找到最大的检验数所在位置disp(n1);disp("--------------");if all(A(:,n1) <= 0)%判断无界解,否则继续迭代x_opt=[];flag=3;break;endb1 = 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;  %基矩阵的逆乘以bA(:,bs) = I;  %基矩阵更新为单位矩阵
end
if flag==1disp('唯一最优解');return
elseif flag==2disp('无穷多最优解');return
elseif flag==3disp('有无界解');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,:)])==Iind_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);                %计算cBSigma = zeros(1,n); %检验数数组for i=1:nz=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);returnenditer=iter+1;% 换基ind_B(ind_B == r) = s;      %新的基变量索引ind_N = setdiff(1:n, ind_B); %非基变量索引% 更新A和bA(:,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)

运行结果