实验一:线性规划单纯形算法
一、实验目的
通过实验熟悉单纯形法的原理,掌握Matlab 循环语句的应用,提高编
程的能力和技巧。
二、算法
对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始
基本可行解。设初始基为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 , B wB C =,得到1B w C B -=,对于非基变量,计算判别数
1i i i B i i z c c B p c σ-=-=-,令 max{}k i i i R
z c σ∈=-,R 为非基变量集合 若判别数0k σ≤ ,则得到一个最优基本可行解,运算结束;否则,转到下一
步
(3).解k k By p =,得到1k
k y B p -=;若0k y ≤,即k y 的每个分量均非正数,则停止计算,问题不存在有限最优解,否则,进行步骤(4).
(4).确定下标r,使{}:0min ,0t r
rk tk tk 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 -=-
四、计算框图
是
否
是
否
五、计算程序
function [x,f]=zuiyouhua(A,b,c)
初始可行解B 令1,0,B N B B x B b b x f c x -==== 计算单纯形乘子1B w c B -=,计算判别数,i j j wp c j R σ=-∈(非基变量)令max{,}k j j R σσ=∈ 0?k σ≤ 得到最优解 解方程k k By p =,得到1k k y B p -=。 0?k y ≤ 不存在有限最优解 确定下标r ,是 {}:0min ,0t r rk tk tk b b tk y y t y y >=>且 k x 为进基变量,用k p 替换r B p ,得到新的基矩阵B
size(A)=[m,n];
i=n+1:n+m;%基变量集合,后面m个松弛变量为初始基变量; N=1:n;%初始非基变量;
B=eye(m,m);
xb=b';
xn=zeros(m,1);
f1=0;
w=zeros(1,m);
z=-c;%初始判别数;
flag=1;
while(1)
[a,k]=max(z);%x(k)为进基变量;
if a<=0
flag=0;
break
else
y=inv(B)*A(:,k)
if y<=0
flag=0;
fprintf('不存在最优解')
break
end
t=find(y>0);
[a,r1]=min(b1(t)./y(t))
r=t(r1); %基变量中第r 个变量为退基变量;
i(:,r)=k
B(:,r)=A(:,k);%换基,即将原基中第r个变量换成第k个变量;
cb=c(:,i);%新的价值系数;
xb=inv(B)*b;
b0=xb;
x=zeros(1,n+m)
x(:,i)=xb'
f=cb*xb
z=cb*inv(B)*A-c;%可用z=cb*(B\A)-c,判别数.
end
end
六、数值实验及结果分析
求解线性规划问题:
???????=≥≤-=+-=++--4
,3,2,1,012
216443033..3min 21421
3212
1i x x x x x x x x x t s x x i
在工作区输入:
A=[3,3,1,0;-4,-4,0,1;2,-1,0,0];
b=[30,16,12]';
c=[-3,1,0,0];
[x,f]=zuiyouhua(A,b,c)
x =
7.3333 2.6667 0 0 0 56.0000 0
f =
-19.3333
检验结果正确
实验一:线性规划单纯形算法 一、实验目的 通过实验熟悉单纯形法的原理,掌握Matlab 循环语句的应用,提高编程的能力和技巧。 二、实验用仪器设备、器材或软件环境 Windows Xp 操作系统 ,Matlab6.5,计算机 三、算法 对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始 基本可行解。设初始基为B,然后执行如下步骤: (1).解B Bx b =,求得1 B x B b -=,0,N B B x f c x ==令计算目标函数值 1(1,2,...,)i m B b i -=i 以b 记的第个分量 (2).计算单纯形乘子w , B wB C =,得到1 B w C B -=,对于非基变量,计算判别数 1i i i B i i z c c B p c σ-=-=-,令 max{}k i i i R z c σ∈=-,R 为非基变量集合 若判别数0k σ≤ ,则得到一个最优基本可行解,运算结束;否则,转到下一步 (3).解k k By p =,得到 1 k k y B p -=;若0k y ≤,即k y 的每个分量均非正数,则停止计算,问题不存在有限最优解,否则,进行步骤(4). (4).确定下标r,使 { } :0 min ,0 t r rk tk tk 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,B N B B x B b b x f c x -==== 计算单纯形乘子1 B w c B -=,计算判别数,i j j wp c j R σ=-∈(非基变量) 令max{,}k j j R σσ=∈ 0?k σ≤ 得到最优解 解方程k k By p =,得到1k k y B p -=。 0?k y ≤ 不存在有限最优解 确定下标r ,是 { }:0 min ,0 t r rk tk tk b b tk y y t y y >=>且 k x 为进基变量,用 k p 替换r B p ,得到新的基矩阵B
单纯形法演算 j c 2 1 B C X B b 1x 2x 3x 4x 5x 0 3x 15 0 5 1 0 0 无穷 0 4x 24 6 2 0 1 0 4 0 5x 5 1 1 0 0 1 5 j j z c -(检验数) 2 1 首先列出表格,先确定正检验数最大值所在列为主列,然后用b 除以主列上对应的同行数字。除出来所得值最小的那一行为主行,根据主行和主列可以确定主元(交点)。接着把主元化为1并把X4换成X1. ??? ??? ?≥=++=++=+++++=0,,524261550002max 5152 14213 25 4321x x x x x x x x x x x x x x x z ??????? ≥≤+≤+≤+=0 ,5 24261552max 21212122 1x x x x x x x x x z
j c 2 1 B C X B b 1x 2x 3x 4x 5x 0 3x 15 0 5 1 0 0 2 1x 4 1 2/6 0 1/6 0 0 5x 5 1 1 0 0 1 j j z c - 2 1 这时进行初等行列变换,把主列换单位向量,主元为1。也就是X5所在行减去X1所在行。并且重新计算检验数。 j c 2 1 B C X B b 1x 2x 3x 4x 5x 0 3x 15 0 5 1 0 0 2 1x 4 1 2/6 0 1/6 0 0 5x 5-4 1-1=0 1-2/6 =4/6 0-1/6=-1/6 1 j j z c - 2-2*1-0*0-0*1=0 1-0*5-2*2/6-0*4/6=1/3 0-0*0-2*1/6-0*-1/6=-1/3 再次确定主元。为4/6。然后把X5换成X2。并且把主元化成1。
数 学 软 件 与 实 验 数学与信息科学学院 信息与计算科学
单纯形法的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:n if sgm(i)>0 h=1; end end while h>0 [msg,mk]=max(sgm); for i=1:m sta(i)=b(i)/A(i,mk); end [mst,mr]=min(sta); zy=A(mr,mk); for i=1:m
if i==mr for j=1:n A(i,j)=A(i,j)/zy; end b(i)=b(i)/zy; end end for i=1:m if i~=mr for j=1:n A(i,j)=A(i,j)-A(i,mk)*A(mr,j); end b(i)=b(i)-A(i,mk)*b(mr); end end B1=A(:,1:m); cb(mr)=c(mk); xx(mr)=mk; sgm=c-cb*B1*A; for i=m+1:n if sgm(i)>0 h=1;
end end end fm=c*xx; 例题: 编写下列求解如下线性规划问题的单纯形法函数min f'x s.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:l d=find(c(:,lie(j)));
%求解标准型线性规划: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 flag kk=kk+1; if A(mA,:)<=0 %已找到最优解 flag=0; sol=zeros(1,nA-1);%给每个变量赋初值0 for i=1:mA-1 sol(N(i))=A(i,nA);%给基变量赋新值(替换0) end %给出最优解 val=-A(mA,nA); else for i=1:nA-1 if A(mA,i)>0&A(1:mA-1,i)<=0 %问题有无界解 disp('have infinite solution!'); flag=0; break; end end if flag %还不是最优表,进行转轴运算 temp=0; for i=1:nA-1 if A(mA,i)>temp temp=A(mA,i); inb=i; % 进基变量的下标 end end %选择最大检验数纵向对应的变量为进基变量 sita=zeros(1,mA-1); for i=1:mA-1 if A(i,inb)>0 sita(i)=A(i,nA)/A(i,inb); end end temp=inf; for i=1:mA-1 if sita(i)>0&sita(i) 算法实现与分析 算法1.单纯形法 具体算例: 标准化后: 用单纯形法求解,程序如下: clear clc M=1000000; A=[3,2,-3,1,0;1,-2,1,0,1];%系数矩阵 C=[-3,1,2,M,M,0];%价值矩阵 B=[6;4]; Xt=[4 5]; for i=1:length(C)-1 D=0; for j=1:length(Xt) D=D+A(j,i)*C(Xt(j)); end xi(i)=C(i)-D; end s=[]; for i=1:length(xi) if xi(i)<0 s=[s,i]; end end f=length(s); h=1; while(f) for k=1:length(s) j=1; A x=[]; for i=1:length(Xt) if A(i,s(k))>0 x(j)=i; j=j+1; end end x if(length(x)+1==1) break; end y=1 x for i=1:length(x) if B(x(i))/A(x(i),s(k)) 单纯形法(Matlab 程序) %%单纯形法( Matlab 程序) a=input('input the major matrix A '); b=input('input the matrix b '); n=input('input the judgement '); %%为计数器(确定循环次数) g=0; while g<40 %%确定非负 alength=max(size(n)); blength=max(size(b)); m=0; for i=1:alength if n(i)>=0 m=m+1; end end; if m==alength x=b; break end; %%找 K s=min(n); for i=1:alength if n(i)==s k=i; break end; end; %%a[i,k] 的非负性 m=0; for i=1:blength if a(i,k)<0 m=m+1; end; end; if m==blength disp('x does not exit'); judge=1; break end; %%找 L 确定主元 cc=100000; for i=1:blength if a(i,k)>0 if (b(i)/a(i,k))< cc cc=b(i)/a(i,k); end end end; for i=1:blength if a(i,k)~=0 if (b(i)/a(i,k))==cc l=i; break end end end; %%计算 ,a 标准化 zu=a(l,k); aa=a; for i=1:l-1 for j=1:alength aa(i,j)=a(i,j)- a(l,j)*a(i,k)/a(l,k); end end; for i=l+1:blength for j=1:alength aa(i,j)=a(i,j)- a(l,j)*a(i,k)/a(l,k); end end; for j=1:alength aa(l,j)=a(l,j)/zu; end; %%b勺判别 bb=b; bb(l)=b(l)/zu; for i=1:l-1 bb(i)=b(i)- b(l)*a(i,k)/a(l,k); end; for i=l+1:blength bb(i)=b(i)- b(l)*a(i,k)/a(l,k); end; b=bb; %%确定判别数 线性规划 线性规划是处理线性目标函数和线性约束的一种较为成熟的方法,目前已经广泛应用于军事、经济、工业、农业、教育、商业和社会科学等许多方面。 8.2.1 基本数学原理 线性规划问题的标准形式是: ????? ??????≥=+++=+++=++++++=0,,,min 21221122222121112 121112211n m n mn m m n n n n n n x x x b x a x a x a b x a x a x a b x a x a x a x c x c x c z 或 ???? ?????=≥===∑∑==n j x m i b x a x c z j n j i j ij n j j j ,,2,1,0,,2,1,min 1 1 写成矩阵形式为: ?? ???≥==O X b AX CX z min 线性规划的标准形式要求使目标函数最小化,约束条件取等式,变量b 非负。不符合这几个条件的线性模型可以转化成标准形式。 MATLAB 采用投影法求解线性规划问题,该方法是单纯形法的变种。 8.2.2 有关函数介绍 在MATLAB 工具箱中,可用linprog 函数求解线性规划问题。 linprog 函数的调用格式如下: ●x=linprog(f,A,b):求解问题minf'*x ,约束条件为A*x<=b 。 ●x=linprog(f,A,b,Aeq,beq):求解上面的问题,但增加等式约束,即Aeq*x=beq 。若没有不等式约束,则令A=[ ],b=[ ]。 ●x=linprog(f,A,b,Aeq,beq,lb,ub):定义设计x 的下界lb 和上界ub ,使得x 始终在该范围内。若没有等式约束,令Aeq=[ ],beq=[ ]。 ●x=linprog(f,A,b,Aeq,beq,lb,ub,x0):设置初值为x0。该选项只适用于中型问题,默认时大型算法将忽略初值。 ●x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options):用options 指定的优化参数进行最小化。 ●[x,fval]=linprog(…):返回解x 处的目标函数值fval 。 ●[x,lambda,exitflag]=linprog(…):返回exitflag 值,描述函数计算的退出条件。 ●[x,lambda,exitflag,output]=linprog(…):返回包含优化信息的输出参数output 。 ●[x,fval,exitflag,output,lambda]=linprog(…):将解x 处的拉格朗日乘子返回到lambda 参数中。 科学出版社《运筹学》教材 第一章引言 第二章线性规划,姜林 第三章对偶规划,姜林 第四章运输问题,姜林 第五章整数规划,姜林 第六章非线性规划,姜林 第七章动态规划,姜林 第八章多目标规划,姜林 第九章图与网络分析,熊贵武 第十章排队论,熊贵武 第十一章库存论,王勇 第十二章完全信息博弈,王勇 第十三章不完全信息博弈,王勇 第十四章决策论与影响图 第十五章运筹学模型的计算机求解 成年人每天需要从食物中摄取的营养以及四种食品所含营养和价格见下表。问 如何选择食品才能在满足营养的前提下使购买食品的费用最小? 食品名称热量(kcal) 蛋白质(g) 钙(mg)价格(元)猪肉1000 50 400 14 鸡蛋800 60 200 6 大米900 20 300 3 白菜200 10 500 2 营养需求量 2000 55 800 解:设需猪肉、鸡蛋、大米和白菜各需 x1,x2,x3,x4斤。则热量的需求量为: 2000 20090080010004 3 2 1 x x x x 蛋白质 某工厂要做100套钢架,每套有长 3.5米、2.8米和2根2.4米的圆钢组成(如右图)已知原 料长12.3米,问应如何下料使需用的原材料最省。 解:假设从每根 12.3米的原材料上截取 3.5米、2.8米和2根2.4 米,则每根原材料需浪费 1.2米,做100套需浪费材料 120米,现 采用套裁的方法。 方案一二三四五六3.5 2.8 2.4 0 0 5 0 4 0 1 2 1 1 3 0 2 0 2 2 1 1 合计剩余 12 0.3 11.2 1.1 11.5 0.8 11.9 0.4 11.8 0.5 12.2 0.1 现在假设每种方案各下料x i (i=1、2、3、4、5、6),则可列出方程: minZ=0.3x 1+1.1x 2+0.8x 3+0.4x 4+0.5x 5+0.1x 6 约束条件: x 3+x 4+2x 5+2x 6=100 4x 2+2x 3+3x 4+x 6=100 5x 1+x 3+2x 5+x 6=200 ,,,800 50030020040055 102060503000 2009008001000. .23614min 4 3214 3 2 1 4 32 14 32 14321x x x x x x x x x x x x x x x x t s x x x x z 多目标线性规划的若干解法及MATLAB 实现 一.多目标线性规划模型 多目标线性规划有着两个和两个以上的目标函数,且目标函数和约束条件全是线性函 数,其数学模型表示为: 11111221221122221122max n n n n r r r rn n z c x c x c x z c x c x c x z c x c x c x =+++??=+++?? ??=+++? (1) 约束条件为: 1111221121122222112212,,,0 n n n n m m mn n m n a x a x a x b a x a x a x b a x a x a x b x x x +++≤??+++≤?? ??+++≤?≥?? (2) 若(1)式中只有一个1122i i i in n z c x c x c x =+++ ,则该问题为典型的单目标线性规划。我们记:()ij m n A a ?=,()ij r n C c ?=,12(,,,)T m b b b b = ,12(,,,)T n x x x x = , 12(,,,)T r Z Z Z Z = . 则上述多目标线性规划可用矩阵形式表示为: max Z Cx = 约束条件:0 Ax b x ≤?? ≥? (3) 二.MATLAB 优化工具箱常用函数[3] 在MA TLAB 软件中,有几个专门求解最优化问题的函数,如求线性规划问题的linprog 、求有约束非线性函数的fmincon 、求最大最小化问题的fminimax 、求多目标达到问题的fgoalattain 等,它们的调用形式分别为: ①.[x,fval]=linprog(f,A,b,Aeq,beq,lb,ub) f 为目标函数系数,A,b 为不等式约束的系数, Aeq,beq 为等式约束系数, lb,ub 为x 的下 限和上限, fval 求解的x 所对应的值。 算法原理:单纯形法的改进方法投影法 ②.[x,fval ]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub ) fun 为目标函数的M 函数, x0为初值,A,b 为不等式约束的系数, Aeq,beq 为等式约束 北京联合大学 实验报告 项目名称:运筹学专题实验报告 学院:自动化专业:物流工程 班级: 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 cx 1.求解上述一般标准线性规划的单纯形算法(修正)步骤如下: 对于一般的标准形式线性规划问题(求极小问题),首先给定一个初始基本可行解。设初始基为B,然后执行如下步骤: (1).解B Bx b =,求得 1 B x B b -=,0,N B B x f c x ==令计算目标函数值 1(1,2,...,)i m B b i -=i 以b 记的第个分量 (2).计算单纯形乘子w, B wB C =,得到1 B w C B -=,对于非基变量,计算判别 数1i i i B i i z c c B p c σ-=-=-,可直接计算 σ =1 B A c c B --令 max{}k i R σσ∈=,R 为非基变量集合 若判别数0k σ≤ ,则得到一个最优基本可行解,运算结束;否则,转到下一 步 (3).解k k By p =,得到 1 k k y B p -=;若0k y ≤,即k y 的每个分量均非正数, 则停止计算,问题不存在有限最优解,否则,进行步骤(4).确定下标r,使 { }:0 min ,0 t r rk tk tk b b tk y y t y y >=>且r B x 为离基变量, ,r k B x p k 为进基变量,用p 替换得到新的基矩阵B,还回步骤(1) ; 大连民族学院 数学实验报告 课程:____________________ 最优化方法______________________ 实验题目: ___________ 单纯形法的matlab实现 ___________________ 系别:______________________ 理学院________________________ 专业:__________________ 信息与计算科学____________________ 姓名:__________________________________________________ 班级:_____________________ 信息102班 ____________________ 指导教师:___________________ 葛仁东_______________________ 完成学期:2013 年__9 ___________ 月_2________ 日 实验目的: 实验方法和步骤(包括数值公式、算法步骤、程序) : 考察标准形式的线性规划问题: min f(x) C T x s.t Ax b, x 0 设x(k)F为一个基本可行解,单纯形方法首先检验它的最优性。如果它不是最优的,确定与该顶点相连的一条使目标函数下降的边;接下来确定沿这个边移 动多远可以到达另一个更优的相邻点,也就是得出一个新的基本可行解 算法步骤: 步骤1给定一个初始基本可行解,记迭代次数 k 1 ; 步骤2 :计算单纯形乘子y k B k T c B k)和简约价值系数向量C N k) c N k) N T y k ; 步骤3 :最优性检验,计算C?k) min{C (k)|j 2},如果C?k) 0,则x (k)为最优解, 停止迭代;否则有x p 0,选x p 为入基变量; 步骤4:确定出基变量,计算g k) B k 1a p ,如果对所有j B k ,有器)0,则问题 无有界的最优解,停止迭代;否则确定出基变量指标 步骤5:交换B k 的列a q 与N k 的列a p 得到新的基矩阵 盼和山+1,计算新的基本可 行解 x (k1),置k:k 1后转步骤 2; 在上述算法中,当存在不止一个简约价值系数 C j k) 0时,选取最负的?“的 指标为p ,并以X p 作为入基变量。 Matlab 计算程序: Function] x,f]=zuiyouhua(A,b,c) Size(A)=[m, n]; i=n+1: n+m; N=1: n; B=eye(m,m); xb=b '; xn=zeros(m,1); f1=0; w=zeros(1,m); z=-c; flag=1; while(1) [a,k]=max(z); If a<=0 flag=0; break else y=i nv(B)*A(:,k) b (k) B k }; min{ _(k )殆 0, j function Z=dcxf(c,A,N) %定义函数名称为dcxf。 l=length(N); CB=c(N(1):N(l)) [m,n]=size(A); b=A(:,n); A=A(:,1:n-1); %参数包括目标函数系数(C),约束条件的系数矩阵(A), %其中A的最后一列为约束条件的右端值b,初始基向量的位置(N)。 sigma=c-CB*A;%计算检验数sigma。 display('初始单纯形表为:');%输出初始的单纯形表table=[nan,nan,nan,c;CB',N',b,A;nan,nan,nan,si gma] opt=1;step=0; while opt step=step+1;%定义循环,直到第“step”步找到最优解(opt=0)。 if sum(sigma>0)==0 %利用检验数判断是否得到最优解,并给出提示。 display('没有得到最优解,继续迭代.'); opt=0; else inb=find(sigma==max(sigma)); %利用单纯形方法找到入基变量的位置 num=length(inb); Inb=inb(num) flag=0; for i=1:m %利用单纯形方法找出出基变量的位置 if A(i,inb)>0 theta(i)=b(i)/A(i,inb); else theta(i)=inf; end end outb=find(theta==min(theta)); num=length(outb); %判断足否出现退化现象,如出现退化,给il{语言提示,并取最后出现的符合出基条件的变量为出基变量。 if num~=1 display('出现退化情况.'); end outb=outb(num); 单纯形法(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)〉=0 m二m+1; end end; if m==alength x=b; break end; %%找K s二min(n); for i=l:alength if n(i) ==s k二i; break end; end; %%a[i,k]的非负性 m=0; for i=l:blength if a(i, k)<0 m二m+1; end; end; if m==blength disp('x does not exit'); judge二1; break end; %%找L确定主元cc=100000; for i=l:blength if a (i, k) >0 if (b(i)/a(i, k))< cc cc=b(i)/a(i, k); end end end; for i=l:blength if a(i, k)~=0 if (b(i)/a(i, k))==cc 1二i; break end end 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:blength for j=l :alength aa(i, j)=a(i, j)- a(l, j)*a(i, k)/a(l, k); end end; for j=l:alength aa(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; %%确定判别数 实验报告 实验题目: 单纯形法的matlab实现学生: 学号: 实验时间: 2013-4-15 一.实验名称: 单纯形法的MATLAB 实现 二.实验目的及要求: 1. 了解单纯形算法的原理及其matlab 实现. 2. 运用MATLAB 编辑单纯形法程序解决线性规划的极小化问题, 求出最优解及目标函数值. 三.实验容: 1. 单纯形方法原理: 单纯形方法的基本思想, 是从一个基本可行解出发, 求一个使目标函数值有所改善的基本可行解; 通过不断改进基本可行解, 力图达到最优基本可行解. 对问题 . 0 ,A s.t. def min ≥=x b x cx f 其中A 是一个m ×n 矩阵, 且秩为m, c 为n 维行向量, x 为n 维列向量, b 为m 维非负列向量. 符号“”表示右端的表达式是左端的定义式, 即目标函数f 的具体形式就是cx . 记 ),...,,(n 21p p p A = 令A =(B,N), B 为基矩阵, N 为非基矩阵, 设 ?? ? ???=0B 1-) 0(b x 是基本可行解, 在) 0(x 处的目标函数值 b c b c c cx f 1 -B 1-N B ) 0(0B 0B ),(=?? ????==, 其中B c 是c 中与基变量对应的分量组成的m 维行向量; N c 是c 中与非基变量对应的分量组成的n-m 维行向量. 现由基本可行解) 0(x 出发求解一个改进的基本可行解. 设?? ? ? ??=N B x x x 是任一可行解, 则由b Ax =得到 N 1-1-B N B B x b x -=, 在点x 处的目标函数值 ∑∈--=?? ? ???==R j j j j f x x c c cx f x )c z (),(0N B N B , 其中R 是非基变量下标集, j j p c 1-B B z =. 2. 单纯形方法计算步骤: 首先给定一个初始基本可行解, 设初始基为B, 然后执行下列主要步骤: (1)解b x B =B , 求得_ 1 b b B x B ==-, 令0=N x , 计算目标函数值B B x c f =. (2)求单纯形乘子w , 解B c wB =, 得到1 -=B c w B . 对于所有非基变量, 计算判别数 j j j j j c c -z -=p w . 令 } c -{z max c -z j j R j k k ∈=. 若0c -z k k ≤, 则对于所有非基变量0c -z j j ≤, 对应基变量的判别数总是为零, 因此停止计算, 现行基本可行解是最优解. 否则, 进行下一步. (3)解k k p By =, 得到k 1 k p B y -=, 若0k ≤y , 即k y 的每个分量均非正数, 则停止计算, 问题不存在有限最优解. 否则进行步骤(4). (4)确定下标r, 使 康福 0031 一.无约束优化方法 无约束优化方法的必要性 一般机械优化设计问题,都是在一定的限制条件下追求某一指标为最小,它们都属于约束优化问题。但是为什么要研究无约束优化问题? (1)有些实际问题,其数学模型本身就是一个无约束优化问题。 (2)通过熟悉它的解法可以为研究约束优化问题打下良好的基础。 (3)约束优化问题的求解可以通过一系列无约束优化方法来达到。所以无约束优 化问题的解法是优化设计方法的基本组成部分,也是优化方法的基础。 (4)对于多维无约束问题来说,古典极值理论中令一阶导数为零,但要求二阶可 微,且要判断海赛矩阵为正定才能求得极小点,这种方法有理论意义,但无实用价值。和一维问题一样,若多元函数F(X)不可微,亦无法求解。但古典极值理论是无约束优化方法发展的基础。 共轭梯度法 目前已研究出很多种无约束优化方法,它们的主要不同点在于构造搜索方向上的差别。 (1)间接法——要使用导数,如梯度法、(阻尼)牛顿法、变尺度法、共轭梯度 法等。 (2)直接法——不使用导数信息,如坐标轮换法、鲍威尔法单纯形法等。 用直接法寻找极小点时,不必求函数的导数,只要计算目标函数值。这类方法较适用于解决变量个数较少的(n ≤20)问题,一般情况下比间接法效率低。 间接法除要计算目标函数值外,还要计算目标函数的梯度,有的还要计算其海赛 矩阵。 搜索方向的构成问题乃是无约束优化方法的关键。 1(0,1,2,) k k k k s k α+=+=L x x 共轭梯度法是沿着共轭方向进行搜索,属于共轭方向法中的一种,该方法中每一个共轭向量都是依赖于迭代点处的负梯度而构造出来。共轭梯度法作为一种实用的迭代法,它主要有下面的优点: (1)算法中,系数矩阵A的作用仅仅是用来由已知向量P 产生向量W=AP ,这不仅 可充分利用A的稀疏性,而且对某些提供矩阵A较为困难而由已知向量P 产生向量W=AP 又十分方便的应用问题是很有益的。 (2)不需要预先估计任何参数就可以计算,这一点不像SOR 等; (3)每次迭代所需的计算,主要是向量之间的运算,便于并行化。 共轭梯度法原理的知识较多,请详见《机械优化设计》第四章的第四、五节。 图1为共轭梯度法的程度框图 图1为共轭梯度法的程度框图 二.设计题目及要求 设计题目 用共轭梯度法求二次函数 2 21212 112(,)242f x x x x x x x =+-- 的极小点及极小值。 设计要求 (1)使用matlab 编写程序,熟练撑握matlab 编程方法。 (2)学习并撑握共轭梯度法的原理、方法及应用,并了解不同无约束优化方法的 区别、优缺点及特殊要求。 (3)编写程序,计算出二次函数的极小点及极小值,并适当选取不同的初始点及 迭代精度精度,分析比较结果。 clear clc X=[1 2 3 4 5]; A=[ 1 2 1 0 0; 4 0 0 1 0; 0 4 0 0 1]; C=[2 3 0 0 0 ]; b=[8;16;12]; t=[3 4 5]; B0=A(:,t); while 1 CB0=C(:,t); XN01=X; for i=1:length(t); for j=1:length(X); if XN01(j)==t(i) XN01(j)=0; end end end j=1; for i=1:length(X); if XN01(i)~=0 XN0(j)=XN01(i); j=j+1; end end for j=1:length(XN0); CN0(j)=C(XN0(j)); end N0=[]; for i=1:length(XN0); N0=[N0,A(:,XN0(i))]; end xiN0=CN0-CB0*B0*N0; j=1; z=[]; for i=1:length(xiN0) if xiN0(i)>0 z(j)=i; j=j+1; end end if length(z)+1==1; break; end n=1; for i=1:length(z) if z(i)>z(n) n=i; end end k=XN0(z(n));%换入变量 B=B0*b; P=B0*A(:,k); j=1; for i=1:length(P) if P(i)>0 x(j)=i; j=j+1; end end y=1; for i=1:length(x) if B(x(y))/P(x(y))>B(x(i))/P(x(i)) y=i; end end y1=x(y); y=t(y1);%换出变量 for i=1:length(t) if t(i)==y m=i; break; end end t(m)=k; P2=B0*A(:,k); q=P2(y1); P2(y1)=-1; P2=-P2./q; E=[1 0 0;0 1 0;0 0 1]; E(:,m)=P2; B0=E*B0; end CB0*B0*b 实验:编制《线性规划》计算程序 一、实验目的: (1)使学生在程序设计方面得到进一步的训练,掌握Matlab (C 或VB)语言进行程序设计中一些常用方法。 (2)使学生对线性规划的单纯形法有更深的理解. 二、实验用仪器设备、器材或软件环境 计算机, Matlab R2009a 三、算法步骤、计算框图、计算程序等 本实验主要编写如下线性规划问题的计算程序: ???≥≥≤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, B wB 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 r rk tk tk b b tk y y t y y >=>且r B x 为离基变量, ,r k B x p k 为进基变量,用p 替换得到新的基矩阵B,还回步骤(1); 单纯形法例题 1、例1、目标函数max z=2 * +3 禹+ 2x2 W 8' 4xi W 16 4x2 W 12 k Ki,财鼻0』 解:首先要将约束条件化为标准形:由此可以看出我们需要加上三个松弛变量, ;汁Hi吃:弋"審得到的标准形式为: max z=2~ +3-+ 0 勺+g +O 5 'xt + 2xj + x] = 8 1 4?i X4 =16 4x;+ 巧=12 11 巾弓^3j 乂4, ^5 $ ? 2 3 0 0 0 C R b *4 0 8 1 2 1 0 0 4 0 16 4 0 0 1 0 - 0 ◎12 0[E(|00 1 3 k - z) 2 3 0 0 0 引」一弋木日lk(i才I) 熙=') (也就是如果与主元素同行,则用现在的值除以主元素即可得到即将要填入的值, 否则,就用现在的值减去与主元素构成矩形的边角上的值的乘积再除以主元素之后 的值。例如:上面的第一行所对应的b值为8-(12*2)/4=2 ,故填入值应该为2。而「 则是由我们根据非基变量的检验数的大小,挑选出最大的那个,作为换入变量,然 后用b的值除以该换入变量所在的列的所有值,得到 约束条件: 由于在检验数中仍然存在大于等于的数,而且, 的坐标中有正分量存在,所以需要继续进行迭代运算。通过观察可以看出主元素为1,换入变量为|勒,换出 由于检验数中存在正数,且P5和P3中有正分量存在,所以需要继续迭代(换入变 此时可以发现检验数中没有大于的数,表明已经得到了最优解,所以最优解是: (4,2,0,0,4 ),故目标函数值z=2*4+2*3=14 2、合理利用线材问题,现在要做100套钢架,每套用长为,,和的钢各一根, 单纯形法的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的值直接解决基变量其他为0 A=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:m C(i+n)=-M; Xb(i)=n+i; end Xb=Xb'; Cb=C(1,n+1:n+m); for i=1:n Xn(i)=i; end Xn=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;单纯形法matlab程序
单纯形法MATLAB程序
运用Matlab进行线性规划求解(实例)
单纯形法典型例题
多目标线性规划的若干解法及MATLAB实现
实验二:MATLAB编程单纯形法求解
单纯形法的matlab实现(20200814192014)
单纯形法matlab程序
单纯形法MATLAB程序
单纯形法的matlab实现(极小化问题)
用MATLAB实现共轭梯度法求解实例
改进单纯形法matlab程序
单纯形法C语言程序
单纯形法例题(20210121173229)
单纯形法的matlab编程