单纯形法以及对偶单纯形法的Matlab实现
最近在上《运筹与优化》这门课，讲到了单纯形法这部分，老师让我们上机用Matlab实现单纯形法，数学公式实现过程看不懂的我就暴力模拟了单纯形法的整个过程，包含无界解等情况。
1. 确定初始基矩阵
2. 更新sigmasigmasigma,确定换入基
3. 更新chi_dachi\_dachi_da，确定换出基
4. 判断此时解的情况(最优解、无界解)，是否结束计算过程
5. 初等行变换，更新矩阵系数
A=[2,-3,2,1,0;1/3,1,5,0,1];
b=[15,20]';
c=[1,2,1,0,0];
m=2;
n=5;
format rat;%结果以分数形式显示
[x_opt,fx_opt,iter]=Simplex_eye(A,b,c,n,m)%x_opt为最求解，fx_opt为最优函数值，iter为迭代次数function [x_opt,fx_opt,iter]=Simplex_eye(A,b,c,n,m)%先选出初始单位阵iter=0;init_eye_pos=nchoosek(1:1:n,m);[row,col]=size(init_eye_pos);for i = 1:rowtmp_init_eye=init_eye_pos(i,:);init_eye=A(:,tmp_init_eye);flag=is_ones(init_eye);if flag==1ji_x=tmp_init_eye;break;endend%循环迭代求解sigma=zeros(n);chi_da=(zeros(m))';CB=(zeros(m))';tmp_A=A;flag=1;while flagiter=iter+1;%更新sigma,确定换入基for pos_c = 1:nsigma(pos_c,1)=c(pos_c);for pos_tmp = 1:msigma(pos_c,1)=sigma(pos_c,1)-CB(pos_tmp,1)*tmp_A(pos_tmp,pos_c);endendmax_sigma=-1;for pos_sigma = 1:nif sigma(pos_sigma,1)>max_sigmamax_sigma_pos=pos_sigma;max_sigma=sigma(pos_sigma,1);endendif max_sigma<=0%此时已经找到最优解x_opt=(zeros(n))';x_opt(ji_x,1)=b(:,1);fx_opt=0;for pos_x = 1:nfx_opt=fx_opt+c(pos_x)*x_opt(pos_x,1);endflag=0;else%更新chi_da，确定换出基for chi_da_pos=1:mchi_da(chi_da_pos,1)=b(chi_da_pos,1)/tmp_A(chi_da_pos,max_sigma_pos);endmin_chi_da=-1;for chi_da_pos=1:mif chi_da(chi_da_pos,1)>0 && tmp_A(chi_da_pos,max_sigma_pos)>0%找出min(chi_da)>0的值且对应左边的值要>0(老师说的)if min_chi_da==-1 || chi_da(chi_da_pos,1)<min_chi_da%找到最小的chi_da以及其位置min_chi_da=chi_da(chi_da_pos,1);min_chi_da_pos=ji_x(chi_da_pos);endendendif min_chi_da==-1disp("无界解");%无界解，结束循环flag=0;elseout_x_pos=min_chi_da_pos%换出基in_x_pos=max_sigma_pos%换入基for in_out = 1:mif ji_x(in_out)==out_x_posji_x(in_out)=in_x_pos;CB(in_out,1)=c(in_x_pos);break;endend%对矩阵执行初等变换，因为单纯形法只有换入基的列变化，所以只需将换入基的列化为换出基的列的形式即可%将换入基对应的点值化为1tmp_beishu=tmp_A(in_out,ji_x(in_out));for change_1 =1:ntmp_A(in_out,change_1)=tmp_A(in_out,change_1)/tmp_beishu;endb(in_out,1)=b(in_out,1)/tmp_beishu;%记得将b列的值也要更新%矩阵内部的初等变换for change_pos =1:mif change_pos~=in_out%此列不是换入基beishu=tmp_A(change_pos,in_x_pos);for tt=1:n%暴力更新这一行的值tmp_A(change_pos,tt)=tmp_A(change_pos,tt)-beishu*tmp_A(in_out,tt);endb(change_pos,1)=b(change_pos,1)-beishu*b(in_out,1);%同步更新b的值endendendendend
end
function [flag]=is_ones(a)%判断矩阵是否为单位阵flag=1;[row_a,col_a]=size(a);tmp_a=zeros(row_a);%对应的列只有一个1,其他全是0，否则flag=flase;for i1 =1:row_afor j1 = 1:col_aif a(i1,j1)==1 && tmp_a(j1)==0tmp_a(j1)=1;elseif a(i1,j1)==0continue;elseflag=0;break;endendendfor i1 =1:row_aif tmp_a(i1)==0flag=0;break;endend
end
对偶单纯形法的Matlab实现就是将单纯形法的代码修改了一下，也是根据其单纯形表暴力模拟了实现过程。
1. 确定初始基矩阵
2. 更新bbb，找出min(b)min(b)min(b)，确定换出基
3. 更新sigmasigmasigma，找出sigma/y(i,k)sigma/y_(i,k)sigma/y(i,k)中最小的正数，确定换入基
4. 判断此时解的情况(最优解、无界解)，是否结束计算过程
5. 初等行变换，更新矩阵系数
过程图解：
A=[-1,-2,-1,1,0;-2,1,-3,0,1];
b=[-3,-4]';
c=[-2,-3,-4,0,0]';
m=2;
n=5;
format rat;%结果以分数形式显示
[x_opt,fx_opt,iter]=DSimplex_eye(A,b,c,n,m)%x_opt为最求解，fx_opt为最优函数值，iter为迭代次数
function [x_opt,fx_opt,iter]=DSimplex_eye(A,b,c,n,m)%先选出初始单位阵iter=0;init_eye_pos=nchoosek(1:1:n,m);[row,col]=size(init_eye_pos);for i = 1:rowtmp_init_eye=init_eye_pos(i,:);init_eye=A(:,tmp_init_eye);flag=is_ones(init_eye);if flag==1ji_x=tmp_init_eye;break;endend%循环迭代求解sigma=zeros(n);CB=(zeros(m))';tmp_A=A;flag=1;while flag%更新b，确定换出基min_b=0;for i = 1:mif b(i)<min_bmin_b=b(i);min_b_pos=ji_x(i);min_b_pos_tmp=i;endendif min_b>=0%此时已经找到最优解x_opt=(zeros(n))';x_opt(ji_x,1)=b(:,1);fx_opt=0;for pos_x = 1:nfx_opt=fx_opt+c(pos_x)*x_opt(pos_x,1);endflag=0;break;else%更新sigma,确定换入基for pos_c = 1:nsigma(pos_c,1)=c(pos_c);for pos_tmp = 1:msigma(pos_c,1)=sigma(pos_c,1)-CB(pos_tmp,1)*tmp_A(pos_tmp,pos_c);endendmin_sigma=-1;%找出sigma/y_(i,k)中最小的正数作为换入基，否则无解for i=1:nif i==min_b_poscontinue;endflag_ji=1;for j=1:mif ji_x(j)==iflag_ji=0;break;endendif flag_ji==1&&tmp_A(min_b_pos_tmp,i)<0if min_sigma==-1||min_sigma>sigma(i,1)/tmp_A(min_b_pos_tmp,i)min_sigma=sigma(i,1)/tmp_A(min_b_pos_tmp,i);min_sigma_pos=i;endendendif min_sigma==-1flag=0;disp("存在无界解！");break;endout_x_pos=min_b_pos%换出基in_x_pos=min_sigma_pos%换入基for in_out = 1:mif ji_x(in_out)==out_x_posji_x(in_out)=in_x_pos;CB(in_out,1)=c(in_x_pos);break;endend%对矩阵执行初等变换，因为单纯形法只有换入基的列变化，所以只需将换入基的列化为换出基的列的形式即可%将换入基对应的点值化为1tmp_beishu=tmp_A(in_out,ji_x(in_out));for change_1 =1:ntmp_A(in_out,change_1)=tmp_A(in_out,change_1)/tmp_beishu;endb(in_out,1)=b(in_out,1)/tmp_beishu;%记得将b列的值也要更新%矩阵内部的初等变换for change_pos =1:mif change_pos~=in_out%此列不是换入基beishu=tmp_A(change_pos,in_x_pos);for tt=1:n%暴力更新这一行的值tmp_A(change_pos,tt)=tmp_A(change_pos,tt)-beishu*tmp_A(in_out,tt);endb(change_pos,1)=b(change_pos,1)-beishu*b(in_out,1);%同步更新b的值endendenditer=iter+1;end
end
function [flag]=is_ones(a)%判断矩阵是否为单位阵flag=1;[row_a,col_a]=size(a);tmp_a=zeros(row_a);%对应的列只有一个1,其他全是0，否则flag=flase;for i1 =1:row_afor j1 = 1:col_aif a(i1,j1)==1 && tmp_a(j1)==0tmp_a(j1)=1;elseif a(i1,j1)==0continue;elseflag=0;break;endendendfor i1 =1:row_aif tmp_a(i1)==0flag=0;break;endend
end
单纯形法以及对偶单纯形法的Matlab实现相关推荐
- 单纯形法，修正单纯形法，对偶单纯形法+解释，很长，不过比较完整，对原理解释的比较多，建议有一定心理准备
@TOC最简单的单纯形法是用来解决形如 min /max cx Ax <= b x>=0 这样的问题的,那么如何解决呢,我觉得这篇知乎上的帖子写的很好,直接搬运过来了: 我们直接给出线性规 ...
- C语言实现单纯形法与对偶单纯形法
C语言实现单纯行法与对偶单纯行法 某次为了完成课程要求所做: 单纯形法代码如下,使用方法修改二位数组A[ip][jp], #include<stdio.h> #define ip 3 #d ...
- 运筹说 第31期 | 对偶理论与灵敏度分析—对偶单纯形法
经过前几期的学习,想必大家已经理解了线性规划问题的对偶问题.相关的重要理论以及影子价格.本期,小编将带大家学习对偶单纯形法的计算步骤. 1.对偶单纯形法 (1)单纯形法回顾 在讲解对偶单纯形法之前我们 ...
- 最优化算法对偶单纯形法的matlab实现（对偶单纯形法看这一篇就够了)
文章目录 前言 一.单纯形法表格 1.1可立即读出最优解和最优值的表格具备的特点 二.对偶单纯形法的步骤(流程图) 三.对偶单纯形法的matlab实现 3.1对偶单纯形法matlab代码 3.2测试例 ...
- 运筹学——matlab实现对偶单纯形法
前言 之前我已经发布了matlab实现单纯形法的相关内容,想了解的小伙伴可以点击这个链接: https://blog.csdn.net/abbcdc/article/details/111324457 ...
- 单纯形法详解及MATLAB实现，对偶单纯形法详解及MATLAB实现
单纯形法详解及MATLAB实现,对偶单纯形法详解及MATLAB实现 单纯形法 我们以这样一个方程组做为例子,来看一下单纯形法是如何解题的 这是一个已经化成标准形式的方程组,x4和x5是我们加入的松弛变 ...
- 单纯形法只有两个约束条件_10分钟掌握对偶单纯形法
只听名字的话会感觉对偶单纯形法和对偶问题关系很大,其实不然(想要了解对偶问题的话可以看我之前的文章).对偶单纯形法在我看来和大M法以及两阶段法很像,都是用来补充纯粹的单纯形法无法解决特殊问题的缺陷.而 ...
- 线性规划总结3——单纯形法和对偶单纯形法
单纯形法 1)利用单纯形法原理, 求解原问题 2)求解过程中始终保持常数列≥0(保持原问题可行),而检验数由有正逐渐变为全部≤0(对偶问题由不可行变为可行) 对偶单纯形法 1)利用对偶原理, 求解原问 ...
- 对偶理论和灵敏度分析（单纯形法矩阵形式、对偶理论及转化、影子价格、机会成本、差额成本、对偶单纯形法、灵敏度分析）
对偶理论和灵敏度分析 文章目录 对偶理论和灵敏度分析 修正单纯形法-矩阵描述和计算 对偶问题的提出 原问题与对偶问题转化 对偶问题最优解的经济意义--影子价格 对偶问题约束项的经济意义--生产产品的机 ...
最新文章
- Httpservlet cannot be resolved to a type的原因与解决方法
- 重装windows2003遇到的老问题：0X0000007B和显卡驱动安装不上。
- 计算机基础知识学前自测,2011计算机二级C语言学前自测题：DOS的基本操作
- 骑士聚会问题（PPT截取）
- 查看从库mysql同步错误日志_常见MySQL同步错误恢复方法
- 【HDU - 4342】History repeat itself（数学）
- 【Object C】从Java 一步步走向Object C
- [样本分析] 海莲花CS样本
- 计算机英语pre,如何区别英语前缀pri，pro，per，pre？
- 杂项-公司：Apple
- 面向对象版学员管理系统（存储数据库）
- 程序员如何提升自身价值
- PACS（CT、CR、DR、MR、DSA、RF医院影像管理系统源码）
- java 出路 xls_java生成xls
- 利用Cadence Allegro强大的功能节省您调丝印的时间
- ubuntu linux定制过程
- C++语言对C语言的扩充
- 使用Verilog语言描述计数器——脉动计数器；脉动计数器具有减法计数功能。采用模块设计和行为级设计方法。
- 【ICON查看器】轻量级图标测试工具ICON STRIKE：快速主屏查看
- Qt如何生成license文件