改进单纯形法matlab程序
单纯形法matlab
数学软件与实验数学与信息科学学院信息与计算科学单纯形法的Matlab程序如下:function [xx,fm]=myprgmh(m,n,A,b,c)B0=A(:,1:m);cb=c(:,1:m);xx=1:n;sgm=c-cb*B0^-1*A;h=-1;sta=ones(m,1);for i=m+1:nif sgm(i)>0h=1;endendwhile h>0[msg,mk]=max(sgm);for i=1:msta(i)=b(i)/A(i,mk);end[mst,mr]=min(sta);zy=A(mr,mk);for i=1:mif i==mrfor j=1:nA(i,j)=A(i,j)/zy;endb(i)=b(i)/zy;endendfor i=1:mif i~=mrfor j=1:nA(i,j)=A(i,j)-A(i,mk)*A(mr,j);endb(i)=b(i)-A(i,mk)*b(mr);endendB1=A(:,1:m);cb(mr)=c(mk);xx(mr)=mk;sgm=c-cb*B1*A;for i=m+1:nif sgm(i)>0h=1;endendendfm=c*xx;例题:编写下列求解如下线性规划问题的单纯形法函数min f'xs.t ax<=b(其中b>=0)函数形式function [x,fval,it,op]=singl(f,a,b) 输出中x为最优解fval为最优值it为迭代次数无最优解op=0有最优解op=1编写程序如下:function [x,fval,it,op]=singl(f,a,b)[m,n]=size(a);c=[a eye(m) b;f' zeros(1,m+1)];fval=0;x=zeros(m+n,1);op=1;it=0;e=zeros(1,m);lie=find(f<0);l=length(lie);while(l>0)for j=1:ld=find(c(:,lie(j)));d_l=length(d);if d_l>0for i=1:mif c(i,lie(j))>0e(i)=c(i,end)/c(i,lie(j));elsee(i)=inf;endend[g,h]=min(e);for w=1:m+1if w==hc(w,:)=c(w,:)/c(h,lie(j));elsec(w,:)=c(w,:)-c(h,:)*c(w,lie(j))/c(h,lie(j));endendit=it+1;elseop=0;endendlie=find(c(end,:)<0);l=length(lie);endfor i=1:(m+n)ix=find(c(:,i));if(length(ix)==1)&(ix<=m)&(c(ix,i)==1) x(i)=c(ix,end)elsex(i)=0endendfval=-c(end,end);。
单纯形法MATLAB程序(线性与非线性规划大作业)
单纯形法程序整理.txt [hB,lB]=size(B);%测算B的行、列数 for j=1:hB if(A(B(j,1),B(j,2))~=1)%如果主元不为1,将对应行除以主元化为1 A(B(j,1),:)=A(B(j,1),:)/A(B(j,1),B(j,2)); end for i=1:hA%将其他元化为零 if(i==B(j,1)) continue; end A(i,:)=A(i,:)-A(B(j,1),:)*A(i,B(j,2)); end end end %————————————————% %% 子程序-选主元 function [r,s]=SubFc_XuanZhuYuan(TableK,h,l) %选主元素 % 选出表矩阵TableK中的主元素 %% 选列主元 t=0; for i=1:l if(TableK(h+1,i)<=t) t=TableK(h+1,i);%记录判别数最小值 s=i;%记录下标 end end if(t==0) error('判别数全为正'); end %% 选行主元 for i=1:h%检查元素是否为负,并赋予t初值 if(TableK(i,s)>0) t=TableK(i,l+1)/TableK(i,s); break; end if(i==h)%若元素全为负,则报错 error('主列元素全为非正,规划问题有无界解!'); end end %k=1; for i=1:h %若未知数出现负值,则是初始迭代标准性表产生的,应更换初始主元; if(TableK(i,l+1)<0) error('b向量出现负值,更换初始主元'); end if(TableK(i,s)<=0)%排除A表负元素 continue; end if(TableK(i,l+1)/TableK(i,s)<=t) t=TableK(i,l+1)/TableK(i,s);%记录判别数最小值 r=i;%记录下标 end end %补充判据 k=1; 第 2 页
实验二:MATLAB编程单纯形法求解
北京联合大学实验报告项目名称:运筹学专题实验报告学院:自动化专业:物流工程班级: 1201B 学号:************* 姓名:管水城成绩:2015 年 5 月 6 日实验二:MATLAB 编程单纯形法求解一、实验目的:(1)使学生在程序设计方面得到进一步的训练;,掌握Matlab (C 或VB)语言进行程序设计中一些常用方法。
(2)使学生对线性规划的单纯形法有更深的理解. 二、实验用仪器设备、器材或软件环境 计算机, Matlab R2006三、算法步骤、计算框图、计算程序等本实验主要编写如下线性规划问题的计算程序:⎩⎨⎧≥≥≤0,0..min b x b Ax t s cx 其中初始可行基为松弛变量对应的列组成. 对于一般标准线性规划问题:⎩⎨⎧≥≥=0,0..min b x b Ax t s cx 1.求解上述一般标准线性规划的单纯形算法(修正)步骤如下:对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。
设初始基为B,然后执行如下步骤: (1).解B Bx b=,求得1B x B b -=,0,N B B x f c x ==令计算目标函数值 1(1,2,...,)i m B b i -=i 以b 记的第个分量(2).计算单纯形乘子w,BwB C =,得到1B wC B -=,对于非基变量,计算判别数1i i i B i i z c c B p c σ-=-=-,可直接计算σ=1B A c c B --令max{}k i Rσσ∈=,R 为非基变量集合若判别数0k σ≤ ,则得到一个最优基本可行解,运算结束;否则,转到下一步 (3).解k kBy p =,得到1k k y B p -=;若0k y ≤,即k y 的每个分量均非正数, 则停止计算,问题不存在有限最优解,否则,进行步骤(4).确定下标r,使{}:0min ,0t rrktktk b b tk y y t y y >=>且rB x 为离基变量,,r k B x p k 为进基变量,用p 替换得到新的基矩阵B,还回步骤(1);2、计算框图为:图1 3.计算程序(Matlab):A=input('A=');b=input('b=');c=input('c=');format rat%可以让结果用分数输出[m,n]=size(A);E=1:m;E=E';F=n-m+1:n;F=F';D=[E,F]; %创建一个一一映射,为了结果能够标准输出X=zeros(1,n); %初始化Xif(n<m) %判断是否为标准型fprintf('不符合要求需引入松弛变量')flag=0;elseflag=1;B=A(:,n-m+1:n); %找基矩阵cB=c(n-m+1:n); %基矩阵对应目标值的cwhile flagw=cB/B; %计算单纯形乘子,cB/B=cB*inv(B),用cB/B的目的是,为了提高运行速度。
实验二:MATLAB编程单纯形法求解
北京联合大学实验报告项目名称:运筹学专题实验报告学院:自动化专业:物流工程班级: 1201B 学号:2012100358081 姓名:管水城成绩:2015 年 5 月 6 日实验二:MATLAB 编程单纯形法求解一、实验目的:(1)使学生在程序设计方面得到进一步的训练;,掌握Matlab (C 或VB)语言进行程序设计中一些常用方法。
(2)使学生对线性规划的单纯形法有更深的理解. 二、实验用仪器设备、器材或软件环境 计算机, Matlab R2006三、算法步骤、计算框图、计算程序等本实验主要编写如下线性规划问题的计算程序:⎩⎨⎧≥≥≤0,0..min b x b Ax t s cx 其中初始可行基为松弛变量对应的列组成. 对于一般标准线性规划问题:⎩⎨⎧≥≥=0,0..min b x b Ax t s cx1.求解上述一般标准线性规划的单纯形算法(修正)步骤如下:对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。
设初始基为B,然后执行如下步骤: (1).解B Bx b=,求得1B x B b -=,0,N B B x f c x ==令计算目标函数值 1(1,2,...,)i m B b i -=i 以b 记的第个分量(2).计算单纯形乘子w,BwB C =,得到1B wC B -=,对于非基变量,计算判别数1i i i B i i z c c B p c σ-=-=-,可直接计算σ=1B A c c B --令max{}k i Rσσ∈=,R 为非基变量集合若判别数0k σ≤ ,则得到一个最优基本可行解,运算结束;否则,转到下一步 (3).解k kBy p =,得到1k k y B p -=;若0k y ≤,即k y 的每个分量均非正数, 则停止计算,问题不存在有限最优解,否则,进行步骤(4).确定下标r,使{}:0min ,0t rrktktk b b tk y y t y y >=>且rB x 为离基变量,,r k B x p k 为进基变量,用p 替换得到新的基矩阵B,还回步骤(1);2、计算框图为:图1 3.计算程序(Matlab):A=input('A=');b=input('b=');c=input('c=');format rat%可以让结果用分数输出[m,n]=size(A);E=1:m;E=E';F=n-m+1:n;F=F';D=[E,F]; %创建一个一一映射,为了结果能够标准输出X=zeros(1,n); %初始化Xif(n<m) %判断是否为标准型fprintf('不符合要求需引入松弛变量')flag=0;elseflag=1;B=A(:,n-m+1:n); %找基矩阵cB=c(n-m+1:n); %基矩阵对应目标值的cwhile flagw=cB/B; %计算单纯形乘子,cB/B=cB*inv(B),用cB/B的目的是,为了提高运行速度。
matlab 单纯形法 -回复
matlab 单纯形法-回复什么是单纯形法?单纯形法是一种用于解决线性规划问题的算法。
它是由美国数学家乔治·达特兰德于1947年提出的。
单纯形法通过在一个多面体内移动一个特殊的高维图形,即单纯形,来搜索最优解。
这个多面体被称为可行解区域,而单纯形则是由多个顶点组成的多面体。
单纯形法的步骤如下:1. 理解线性规划问题:在开始使用单纯形法之前,首先需要明确线性规划问题的定义和约束条件。
线性规划问题是一类在线性目标函数下,通过一组线性不等式和等式约束来找到使目标函数最优化的变量值。
2. 转换为标准形式:将线性规划问题转换为标准形式,即将目标函数和约束条件都转化为等式形式,并引入人工变量和松弛变量。
3. 初始化:寻找一个可行基础解,并计算对应的目标函数值。
如果找不到初始的基础解,则问题无解。
4. 选择进基变量:从非基变量中选择一个进基变量,即将其从0增加为正值,以使目标函数值增加。
5. 选择出基变量:确定一个出基变量,即将其从正值减少到0,以保持其他约束条件不变。
6. 单纯形迭代:通过计算目标函数值和约束条件来确定下一个基础解。
如果满足终止条件,则找到了最优解。
否则,继续选择进基变量和出基变量进行下一次迭代。
7. 终止:根据终止条件判定是否找到了最优解。
终止条件可以是找不到进基变量或者出基变量,或者目标函数值已达最优。
单纯形法的优点在于它是一种相对快速和可靠的算法,可以在多项式时间内找到最优解。
然而,它的缺点是当问题规模较大时,计算复杂度会呈指数增长,并且在某些特殊情况下可能会出现挂起、无法终止的情况。
为了改进单纯形法的性能,研究人员提出了一些变种和改进的算法。
例如,内点法通过在多面体内部搜索最优解,避免了单纯形法中移动单纯形的过程,并取得了较好的效果。
其他算法如双轨法、变尺度法等也都在不同程度上改进了单纯形法的缺点。
总而言之,单纯形法是一种经典且有效的用于解决线性规划问题的算法。
它的基本思想是通过迭代逐步优化目标函数值,直到找到最优解。
单纯形法MATLAB程序
单纯形法(Mat lab程序)%%单纯形法(Mat lab程序)a= input (' input the major matrix A '); b=input (' input the matrix b '); n=input C input the judgement ');%%为计数器(确定循环次数)萨0;while g<40%%确定非负alength=max(size(n));blength二max(size(b));m=0;for i=l:alength辻n(i)〉=0m二m+1;endend;if m==alengthx=b;breakend;%%找Ks二min(n);for i=l:alengthif n(i) ==sk二i;breakend;end;%%a[i,k]的非负性m=0;for i=l:blengthif a(i, k)<0m二m+1;end;end;if m==blengthdisp('x does not exit');judge二1;breakend;%%找L确定主元cc=100000;for i=l:blengthif a (i, k) >0if(b(i)/a(i, k))<cccc=b(i)/a(i, k);endend end; for i=l:blengthif a(i, k)~=0if (b(i)/a(i, k))==cc1二i;breakendend end; %%计算,a 标准化zu=a(l, k); aa=a; for i=l:1-1 for j=l:alength aa(i, j)=a(i, j)-a(l, j)*a(i, k)/a(l, k);end end; for i=l+l:blengthfor j=l :alength aa(i, j)=a(i, j)-a(l, j)*a(i, k)/a(l, k);end end; for j=l:alengthaa(l, j)=a(l, j)/zu; end;%%b 勺判别bb=b; bb(l)=b(l)/zu;for i=l: 1~1 bb(i)=b(i)~b⑴*a(i, k)/a(l, k);end;for i=l+l:blength bb(i)二b(i)-b(l)*a(i, k)/a(l, k);end;b二bb; %%确定判别数tt 二n;for j=l:alength11 (j) =n(j)-a(1, j)*n(k)/a(1, k) ; end; n=tt;a=aa;%%显示单纯形表sa sa二[b' aa;0 n];dispC单纯表示例’);disp(g+1);disp(sa);g二g+l;judge=2;end;if judge==2q二0; result=zeros (alength, 2); for j=l+q:alengthif n(j)=0 t=a(:, j) ; zu=find( t) ; resu lt( j, l)=j ; result (j, 2)=x(zu) ; q 二q+1 ;endif n(j)>0 result(j,l)=q+l; q=q+l;endend;dispC最优解’);disp (result);dispC循环次数');end。
实验二MATLAB编程单纯形法求解
北京联合大学实验报告项目名称: 运筹学专题实验报告学院: 自动化专业: 物流工程班级: 1201B 学号:21姓名: 管水城成绩: 2015 年 5 月 6 日实验二:MATLAB 编程单纯形法求解一、实验目的:(1)使学生在程序设计方面得到进一步的训练;,掌握Matlab (C 或VB)语言进行程序设计中一些常用方法。
(2)使学生对线性规划的单纯形法有更深的理解、二、实验用仪器设备、器材或软件环境计算机, Matlab R2006三、算法步骤、计算框图、计算程序等本实验主要编写如下线性规划问题的计算程序:⎩⎨⎧≥≥≤0,0..min b x b Ax t s cx其中初始可行基为松弛变量对应的列组成、对于一般标准线性规划问题:⎩⎨⎧≥≥=0,0..min b x b Ax t s cx1.求解上述一般标准线性规划的单纯形算法(修正)步骤如下:对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。
设初始基为B,然后执行如下步骤:(1)、解B Bx b =,求得1B x B b -=,0,N B B x f c x ==令计算目标函数值 1(1,2,...,)i m B b i -=i 以b 记的第个分量(2)、计算单纯形乘子w, BwB C =,得到1B w C B -=,对于非基变量,计算判别数1i i i B i i z c c B p c σ-=-=-,可直接计算σ=1B A c c B --令 max{}k i R σσ∈=,R 为非基变量集合若判别数0k σ≤ ,则得到一个最优基本可行解,运算结束;否则,转到下一步(3)、解k k By p =,得到1k k y B p -=;若0k y ≤,即k y 的每个分量均非正数, 则停止计算,问题不存在有限最优解,否则,进行步骤(4)、确定下标r,使 {}:0min ,0t rrk tk tk b b tk y y t y y >=>且r B x 为离基变量,,r k B x p k 为进基变量,用p 替换得到新的基矩阵B,还回步骤(1);2、计算框图为:图1 3.计算程序(Matlab):A=input('A=');b=input('b=');c=input('c=');format rat%可以让结果用分数输出[m,n]=size(A);E=1:m;E=E';F=n-m+1:n;F=F';D=[E,F]; %创建一个一一映射,为了结果能够标准输出X=zeros(1,n); %初始化Xif(n<m) %判断就是否为标准型fprintf('不符合要求需引入松弛变量')flag=0;elseflag=1;B=A(:,n-m+1:n); %找基矩阵cB=c(n-m+1:n); %基矩阵对应目标值的cwhile flagw=cB/B; %计算单纯形乘子,cB/B=cB*inv(B),用cB/B的目的就是,为了提高运行速度。
优化设计-单纯形法
x1[0]=0.5*(x1[0]+xe[0]);
x1[1]=0.5*(x1[1]+xe[1]);
x2[0]=0.5*(x2[0]+xe[0]);
x2[1]=0.5*(x2[1]+xe[1]);
}
float max(float x,float y)
{if (x>y)return x;
else {fh=f0; fe=f1; xh=x0; xe=x1; p[0]=0;p[1]=1;p[2]=1;} }
else { fh=f2; fe=f1; xh=x2; xe=x1; p[0]=1;p[1]=1;p[2]=0;} }
else { if(f0<f2)
{ if(f1<f2){fh=f2; fe=f0; xh=x2; xe=x0; p[0]=1;p[1]=1;p[2]=0;}
step3:
{ if(fn4>fh) {sx(); goto step4;}
else {xh[0]=xn4[0];xh[1]=xn4[1];fh=f(xh); goto step4;} }
step4:
eh();
printf("%d\t[%.4f%.4f][%.4f%.4f] [%.4f %.4f] %.4f\n",k,x0[0],x0[1],x1[0],x1[1],x2[0],x2[1],fn1);
1)用单纯形法求法min(x12+2x22-4x1-2x1x2),已知α=1,β=0.5,γ=2,ε=0.005。
程序如下:
#include "stdio.h"
#include "math.h"
实验二:MATLAB编程单纯形法求解
实验二:MATLAB编程单纯形法求解北京联合大学实验报告项目名称:运筹学专题实验报告学院:自动化专业:物流工程班级: 1201B 学号:2021100358081 姓名:管水城成绩:2021 年 5 月 6 日实验二:MATLAB编程单纯形法求解一、实验目的:(1)使学生在程序设计方面得到进一步的训练;,掌握Matlab (C或VB)语言进行程序设计中一些常用方法。
(2)使学生对线性规划的单纯形法有更深的理解. 二、实验用仪器设备、器材或软件环境计算机, Matlab R2021三、算法步骤、计算框图、计算程序等本实验主要编写如下线性规划问题的计算程序:mincx?Ax?b s.t.??x?0,b?0其中初始可行基为松弛变量对应的列组成. 对于一般标准线性规划问题:mincx?Ax?b s.t.??x?0,b?01.求解上述一般标准线性规划的单纯形算法(修正)步骤如下:对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。
设初始基为B,然后执行如下步骤:?1Bx?bx?Bb,令xN?0,计算目标函数值f?cBxB BB(1).解,求得以bi(i?1,2,...,m)记B?1b的第i个分量?1wB?Cw?CBBB(2).计算单纯形乘子w, ,得到,对于非基变量,计算判别数?i?zi?ci?cBB?1pi?ci,可直接计算i?R??cBBA?c令?1 ?k?max{?},R为非基变量集合若判别数步?1By?py?Bpk;若yk?0,即yk的每个分量均非正数, kkk(3).解,得到?k?0 ,则得到一个最优基本可行解,运算结束;否则,转到下一则停止计算,问题不存在有限最优解,否则,进行步骤(4).确定下标r,使bryrk?min?t:ytk?0btytk,且ytk?0?xB为离基变量,rxk为进基变量,用pk替换pBr,得到新的基矩阵B,还回步骤(1);2、计算框图为:开始初始可行基B x?B?1b,x?0,f?cx 令BNBB计算单纯性乘子w?cBB?1,计算判别数?j?wpj?cj,j?R(非基变量)令?k?max{?j,j?R}?k?0是 ? 否得到最优解方程Byk?pk,得到yk?B?1pk,是 yk?0? 否不存在有限确定下标r,使得 ?b?br?min?i|yik?0? yrk?yik?xBr为退基变量,xk进基变量,以pk代替pBr,得到新的基矩阵B 图13.计算程序(Matlab):A=input('A=');b=input('b='); c=input('c=');format rat %可以让结果用分数输出[m,n]=size(A);E=1:m;E=E'; F=n-m+1:n;F=F';D=[E,F]; %创建一个一一映射,为了结果能够标准输出X=zeros(1,n); %初始化X感谢您的阅读,祝您生活愉快。
改进单纯形法寻优的MATLAB实现(1)
6 改进单纯形法的寻优原理
676 改进单纯形法简介
单 纯 形 法 是 应 用 规 则 的 几 何 图 形 !通 过 计 算 单 纯 形 顶 点 的 函 数 值 !根 据 函 数 值 大 小 的 分 布 来 判 断 函 数 变 化 的 趋 势 !然 后 按 一 定 的 规 则 搜 索 寻 优 的 方 法 %8!9’"该 方 法 因 步 长 固 定 !具 有 不 能 加 速 的 缺 点 "改 进 单 纯 形 法 是 在 单 纯 形 法 的 基 础 上 对 步 长 作 适 当 修 改 得 到 的 寻 优 方 法 !在 化 学 化 工 中 应 用 较 广"设需要寻优的目标函数为 :; :<=&!=4!>!=?@!其中 =A<A; &!4!>!?@是自变量!:为响应值"
?@
320175BC
? A /0/2# /0/.1
/
#
D#$/0/#1A/.0/56BC
D$/0//1 .
单纯形法的matlab编程
单纯形法的matlab实现首先输入三个值系数矩阵A目标函数系数行向量C列向量b根据大M法进行扩列A,C,b.使得行数不变,列数增加M 进行的到基向量的坐标,非基变量的坐Cb,Cn,Xb,Xn,此时的值便是典式,不在需要进行进一步化简,只需求解检验变量delta的值迭代过程输入上一步得到A,C,b,Cb,Cn,Xb,Xn,输出值为最优解为X,得到目标函数的最优解Z的值迭代循化用while循环当找到解时结束循环break或者当发现循化结果没有最优解时跳出循环,这里涉及两个判断,两个判断量初始值都可以写在循环外,两者的值共同决定循环的执行与否循化最开始进行判断初始可行解是否为最最优解,若是直接跳出循化,若上面的判断不成立,接下来进行下一个判断,若不符合进行下面入基和出基变量的选值入基和出基变量的循化是两次循化,第一次找到k的值,第二次根据上一次的k找r的值注意因为值有约束,而且是找函数最小值,需要对这个列向量进行变换一下将小于等于0的都变成无穷大,接下来进形下一次的循化,进而找到转轴元将A,b,delta合成一个新的矩阵,进行旋转变化,得到值后反变回相应的值,接下来需要对Xb,Xn的值进行交换这个步骤要两个循环,第一个循化对Ark的所在行进行变化,接下来进行对整个矩阵进行行变换,包括两种情况,两次循化嵌套分别是r==1时和r~=1的时候建立总体X的坐标列向量发生交换时出基变量找Xb,入基变量从X中找有先后顺序先解决Xn的变化。
在解决Xb的值直接解决基变量其他为0A=input('输入系数矩阵\n');b=input('输入列向量b\n');C=input('输入目标函数行向量\n');M=5200;global m;global n;global X;[m,n]=size(A);I=eye(m);A=[A,I];Xb=[];Xn=[];for i=1:mC(i+n)=-M;Xb(i)=n+i;endXb=Xb';Cb=C(1,n+1:n+m);for i=1:nXn(i)=i;endXn=Xn';X=[Xn;Xb];[m,n]=size(A);diedai(A,C,b,Cb,Xb);function[Z]=diedai(A,C,b,Cb,Xb)delta=C-Cb*A;global m;global n;global X;while1s2=0;s1=0;for j=1:nif delta(j)>0s1=1;for i=1:mif A(i,j)>0s2=1;endendendendif s1==0disp('目标函数最优解')Z=Cb*b;disp(Z)disp('基变量为');[Xb,index]=sort(Xb);disp(Xb)b=b(index);disp('基可行解为');disp(b)break;endif s2==0disp('目标函数无界,无最优解');break;end[~,k]=max(delta);p=A(:,k);zhuan=[];for i=1:mzhuan(i)=b(i)/p(i);if zhuan(i)<=0zhuan(i)=inf;endend[~,r]=min(zhuan);b(m+1)=0;Z=[A;delta];Z=[Z,b];z=Z;ark=A(r,k);for j=1:n+1Z(r,j)=Z(r,j)/ark;endif r==1for i=2:m+1for j=1:n+1Z(i,j)=Z(i,j)-z(i,k)*Z(r,j);endendelse for i=[1:r-1,r+1:m+1]for j=1:n+1Z(i,j)=Z(i,j)-z(i,k)*Z(r,j);endendendA=Z(1:m,1:n);delta=Z(m+1,1:n);b=Z(1:m,n+1);Cb(r)=C(k);Xb(r)=X(k);endend。
大连理工大学庞丽萍最优化方法MATLAB程序
班级:优化1班授课老师:庞丽萍姓名:学号:第二章12.(1)用修正单纯形法求解下列LP问题:>>clear>>A=[121100;123010;215001];[m,n]=size(A);b=[10;15;20];r=[-1-2-31];c=[-1-2-31];bs=[3:3];nbs=[1:4];a1=A(:,3);T=A(:,bs);a2=inv(T)*a1;b=inv(T)*b;A=[eye(m),a2];B=eye(m);xb=B\b;cb=c(bs);cn=c(nbs);con=1;M=zeros(1);while conM=M+1;t=cb/B;r=c-t*A;if all(r>=0)x(bs)=xb;x(nbs)=0;fx=cb*xb;disp(['当前解是最优解,minz=',num2str(fx)])disp('对应的最优解为,x=')disp(x)breakendrnbs=r(nbs);kk=find(rnbs==min(rnbs));k=kk(1);Anbs=A(:,nbs);yik=B\Anbs(:,k);xb=B\b;%yi0if all(yik<=0)disp('此LP问题无有限的最优解,计算结束',x)disp(xb)breakelsei=find(yik>0);w=abs(xb(i,1)./yik(i,1));l=find(w==min(w));rr=min(l);yrrk=yik(rr,1);Abs=A(:,bs);D=Anbs(:,k);Anbs(:,k)=Abs(:,rr);Abs(:,rr)=D;F=bs(rr);bs(rr)=nbs(k);nbs(k)=F;AA=[Anbs,Abs];EE=eye(m);EE(:,rr)=-yik./yrrk;Errk=EE;Errk(rr,rr)=1/yrrk;BB=Errk/B;B=inv(BB);cb=c(:,bs);xb=Errk*xb;x(bs)=xb;x(nbs)=0;fx=cb*xb;endif M>=1000disp('此问题无有限最优解')breakendend%结果当前解是最优解,minz=-15对应的最优解为,x=2.5000 2.5000 2.50000第三章30题DFP算法求函数极小点的计算程序function[x,val,k]=dfp(fun,gfun,x0)%功能:用DFP算法求解无约束问题:minf(x)%输入:x0是初始点,fun,gfun分别是目标函数及其梯度%输出:x,val分别是近似最优点和最优值,k是迭代次数.maxk=1e5;%给出最大迭代次数rho=0.55;sigma=0.4;epsilon=1e-5;k=0;n=length(x0);Hk=inv(feval('Hess',x0));%Hk=eye(n);while(k<maxk)gk=feval(gfun,x0);%计算梯度if(norm(gk)<epsilon),break;end%检验终止准则dk=-Hk*gk;%解方程组,计算搜索方向m=0;mk=0;while(m<20)%用Armijo搜索求步长if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk’*dk)mk=m;break;endm=m+1;end%DFP校正x=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x)-gk;if(sk'*yk>0)Hk=Hk-(Hk*yk*yk'*Hk)/(yk'*Hk*yk)+(sk*sk')/(sk'*yk);endk=k+1;x0=x;endval=feval(fun,x0);%习题26的程序调用方式及结果:function y=fun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=(x(1)-1)^2+5*(x2-x(1)^2)^2endfunction y=gfun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[20]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=1.000001.00000val=k=6%习题27的程序调用方式及结果:function y=fun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=x1+2*x(2)^2+exp(x(1)^2+x(2)^2)endfunction y=gfun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[10]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=-0.419360val=0.77291k=536题编写Hooke-Jeeves方法求函数极小点的计算程序。
最优化实验报告(单纯形法的matlab程序,lingo程序)
最优化实验报告(单纯形法的matlab程序,lingo程序)实验一:线性规划单纯形算法一、实验目的通过实验熟悉单纯形法的原理,掌握Matlab 循环语句的应用,提高编程的能力和技巧。
二、实验用仪器设备、器材或软件环境Windows Xp 操作系统 ,Matlab6.5,计算机三、算法对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。
设初始基为B,然后执行如下步骤:(1).解B Bx b =,求得1Bx B b -=,0,N B B x f c x ==令计算目标函数值 1(1,2,...,)i m B bi -=i 以b 记的第个分量(2).计算单纯形乘子w, B wB C =,得到1B wC B -=,对于非基变量,计算判别数1i i i B i i z c c B p c σ-=-=-,令 max{}k i i i Rz c σ∈=-,R 为非基变量集合若判别数0k σ≤ ,则得到一个最优基本可行解,运算结束;否则,转到下一步(3).解k k By p =,得到1k k y B p -=;若0k y ≤,即k y 的每个分量均非正数,则停止计算,问题不存在有限最优解,否则,进行步骤(4). (4).确定下标r,使{}min ,0t rrktktk b b tk y y t y y >=>且r B x 为离基变量。
k x 为进基变量,用k p 替换r B p ,得到新的基矩阵B ,返回步骤(1)。
对于极大化问题,可以给出完全类似的步骤,只是确定进基变量的准则不同。
对于极大化问题,应令min{}k k j j z c z c -=-四、计算框图是否是否开始初始可行解B令1,0,BN B B x B b b x f c x -====计算单纯形乘子1B w c B -=,计算判别数,i j j wp c j R σ=-∈(非基变量)令max{,}kj j R σσ=∈0?k σ≤得到最优解解方程kk By p =,得到1k k y B p -=。
单纯形法的MATLAB实现
clearclcM=1000000;A=[3,2,-3,1,0;1,-2,1,0,1];%约束矩阵C=[-3,1,2,M,M,0];%价值矩阵B=[6,4]';%右端向量s=find(C<0);f=length(s);while(f)for k=1:length(s)x=find(A(:,s(k))>0);y=find(B(x)./A(x,s(1))==min(B(x)./A(x,s(1))));%选择的要有正元素if(length(x)+1==1)break;endendy=x(y);%找到的xj的行数aa=A(y,s(k));%找到的xjA(y,:)=A(y,:)./aa;B(y,:)=B(y,:)./aa;z=find(A(:,s(k)));%除去找到的行z(find(z==y))=[];for i=1:length(z);yz=-A(z(i),s(k));A(z(i),:)=A(z(i),:)+A(y,:)*yz;disp('*')B(z(i),:)=B(z(i),:)+B(y,:).*yz;enddisp('转换后')A=AB=BAB=[A,B];C=C+AB(y,:)*(-C(s(k)))s=find(C<0);vpa([A,B;C]);s=find(C<0);f=length(s);end-C(length(C))%最有解:max 2*x1+3*x2s.t. x1+2*x2<=84*x1<=164*x2<=12x1,x2>=0加入松驰变量,化为标准型,得到A=[1 2 1 0 0 8;4 0 0 1 0 16;0 4 0 0 1 12;2 3 0 0 0 0];N=[3 4 5];然后执行? [sol,val,kk]=ssimplex(A,N)就可以了。
注:基变量对应的基矩阵一定是单位阵。
(这一局限将在后面的升级是改善)% 求解标准型线性规划:max c*x;s.t. A*x=b;x>=0% 本函数中的A是单纯初始表,包括:最后一行是初始的检验数,最后一列是资源向量b % N是初始的基变量的下标%输出变量sol是最优解%输出变量val是最优值,kk是迭代次数function [sol,val,kk]=ssimplex(A,N)[mA,nA]=size(A);kk=0; %迭代次数flag=1;while flagkk=kk+1;if A(mA,:)<=0 % 已找到最优解flag=0;sol=zeros(1,nA-1);for i=1:mA-1sol(N(i))=A(i,nA);endval=-A(mA,nA);elsefor i=1:nA-1if A(mA,i)>0&A(1:mA-1,i)<=0 %? 问题有无界解disp('have infinite solution!');flag=0;break;endendif flag % 还不是最优表,进行转轴运算temp=0;for i=1:nA-1if A(mA,i)>temptemp=A(mA,i);inb=i; % 进基变量的下标endendsita=zeros(1,mA-1);for i=1:mA-1if A(i,inb)>0sita(i)=A(i,nA)/A(i,inb);endendtemp=inf;for i=1:mA-1if sita(i)>0&sita(i)<temptemp=sita(i);outb=i; %出基变量下标endend%以下更新Nfor i=1:mA-1if i==outbN(i)=inb;endend% 以下进行转轴运算A(outb,:)=A(outb,:)/A(outb,inb);for i=1:mAif i~=outbA(i,:)=A(i,:)-A(outb,:)*A(i,inb);endendendendend。