单纯形法详解及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)