最优化实验报告(共20页).doc
精选优质文档-倾情为你奉上数值最优化算法与理论课程实验报告 姓 名： 周飞飞（1） 张琳婧（8） 班 级： 信息与计算科学2班 指导教师：刘陶文2014年11月27日专心-专注-专业数值最优化算法与理论课程实验报告课程名称数值最优化算法与理论班级信息与计算科学2班小组成员周飞飞（1）张琳婧（8）实验课题拟Newton法（BFGS算法）及FR共轭梯度法求解无约束问题实验目的通过上机实验掌握最优化的实用算法的结构及性能，并用这些算法解决实际的最优化问题，掌握一些实用的编程技巧。实验要求选用你喜欢的无约束优化的某种梯度法 (最速下降法，Newton法，拟牛顿法，共轭梯度法)通过编程，上机实验对所提供的测试问题进行测试、运行，然后提供实验报告。在实验报告中指出你选用的算法、参数设置、终止准则、线性搜索以及实验结果，附加你的实验心得。实验内容使用非精确Wolf-Powell线性搜索实现拟牛顿法（BFGS算法）及FR共轭梯度法求解无约束问题，并通过Matlab软件实现算法，观察分析实验过程，对比实验结果来进一步理解两种方法的原理及优点与缺陷。目 录1、 1481213一、 实验原理无约束问题下降算法是求解无约束优化问题的一类最基本的算法。其一般步骤为：（已知近似最优解）1. 首先，计算下降方向满足：2. 然后计算步长满足：3. 计算新的近似最优解：这次实验所运用的拟牛顿法及FR共轭梯度法主要是在下降算法的基础上，求解下降方向的方法上有所不同。（一） 拟牛顿法（1）拟牛顿法的简述 拟(Quasi-Newton Methods)是求解非线性优化问题最有效的方法之一，于20世纪50年代由美国Argonne的物理学家W. C. Davidon所提出来。Davidon设计的这种算法在当时看来是非线性优化领域最具创造性的发明之一。不久R. Fletcher和M. J. D. Powell证实了这种新的算法远比其他方法快速和可靠，使得非线性优化这门学科在突飞猛进。拟牛顿法和(Steepest Descent Methods)一样只要求迭代时知道的梯度。通过测量梯度的变化，构造一个目标函数的模型使之足以产生超线性收敛性。这类方法大大优于最速下降法，尤其对于困难的问题。另外，因为拟牛顿法不需要的信息，所以有时比牛顿法(Newton's Method)更为有效。如今，优化软件中包含了大量的拟牛顿算法用来解决无约束，约束，和大规模的优化问题。（2） 拟牛顿法的思想拟牛顿法是对牛顿法的一种改善保持其优点：快速收敛性克服其缺陷：需计算Hessian矩阵且要求其正定方法：构造对称正定矩阵或使其满足 计算搜索方向：构造的要求如下：（3） 拟牛顿法的BFGS算法BFGS修正公式：该公式由：Broyden-Fletcher-Goldfarb-Shanno 提出，是迄今为止最好的拟牛顿修正公式。BFGS算法：（二） 非线性共轭梯度法（FR算法）（1）共轭梯度法的简述 共轭梯度法是介于与之间的一个方法，它仅需利用信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hessian矩阵并求逆的缺点，不仅是解决大型线性方程组最有用的方法之一，也是解大型最有效的算法之一。 共轭梯度法最早是由Hestenes和Stiefle（1952）提出来的，用于解正定的线性方程组，在这个基础上，Fletcher和Reeves（1964）首先提出了解非线性最优化问题的共轭梯度法。由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现在共轭梯度法已经广泛地应用与实际问题中。 共轭梯度法是一个典型的法，它的每一个搜索方向是互相共轭的，而这些搜索方向d仅仅是负梯度方向与上一次的搜索方向的组合，因此，存储量少，计算方便 。（2） 共轭梯度法的思想Fletcher-Reeves共轭梯度法，简称FR法。共轭梯度法的基本思想是把共轭性与最速下降方法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜素，求出目标函数的极小点。根据共轭方向基本性质，这种方法具有二次终止性。给定初始点，取初始搜索方向，在后面的迭代中取负梯度方向和前一搜索方向的线性组合作为搜索方向即，其中是待定参数，适当选取，使得。FR共轭梯度法，搜索方向的计算公式为： （3） FR共轭梯度法算法由Fletcher-Reeves(1964)提出。FR算法：（三） Wolfe-Powell线搜索Armijo型线性搜索的条件是Wolfe-Powell型线性搜索中的第一个条件。Wolfe-Powell型线性搜索中的第二个条件的作用在于限制过小的步长，减少了求解时的计算量。同时运用非精确Wolfe-Powell线性搜索保证了拟牛顿法（BFDS算法）中的矩阵对称正定。Wolfe-Powell线搜索算法：二、 实验内容拟牛顿法程序：for i=1:nn %1-nn函数依次进入运算(1)初值准备nprob=numer(i); n,m,xk,filename=initf(nprob);% 读初始数据xk=factor*xk;bk=eye(n);k=0;tic; %计时开始 fk=objfcn(n,m,xk,nprob);fnum=1;gk=grdfcn(n,m,xk,nprob);gnum=1;delta=norm(gk,2);(2)迭代开始while k<1000 %迭代上限1000 if delta<=teminate % break;Else(3)确定下降方向 dk=-linsolve(bk+muk*eye(n),gk);%求解下降方向 gk1=gk;fk1=fk;gkdk=gk'*dk; if gk'*dk>=-1.0e-14%当dk不是充分下降时采用负梯度为搜索方向 dk=-gk; end(4)确定步长%利用Wolfe-Powell搜索计算步长 alphak,fk,gk,wfnum,wgnum=wolfe2(n,m,xk,dk,fk1,gk1,nprob);%利用Wolfe-Powell搜索计算步长 (5)计算 fnum=fnum+wfnum; gnum=gnum+wgnum; xk1=xk;xk=xk1+alphak*dk; fk=objfcn(n,m,xk,nprob); gk=grdfcn(n,m,xk,nprob); if norm(gk,2)<=teminate k=k+1; break; end(6)由BFGS修正公式得 % Bk update sk=xk-xk1;bks2=sk'*bk*sk;yk=gk-gk1; yksk=yk'*sk; if yksk>0 bks1=bk*sk*sk'*bk; yks=yk*yk'/yksk;bk1=bk; bk=bk1-bk1*sk*sk'*bk1/(sk'*bk1*sk)+yk*yk'/(yk'*sk); end end k=k+1;End(7)无约束问题运算结束后记录所花费时间time=toc;%终止计时if time<=0. t(i,s)=0.0001;else t(i,s)=time;%将每个无约束问题求解时间记录End(8)输出无约束问题的运行结果fprintf('nt%sttt%2dttt%5dtttt%5dttt%5dttt%4fn',filename,n,k,fnum,gnum,time);%结果输出EndFR共轭梯度法for i=1:nn(1)初值准备nprob=numer(i); n,m,xk,filename=initf(nprob);% 读初始数据xk=factor*xk;bk=eye(n);k=0;tic; %计时开始 fk=objfcn(n,m,xk,nprob);fnum=1;gk=grdfcn(n,m,xk,nprob);gnum=1;delta=norm(gk,2);dk=-gk;(2)迭代开始while k<1000 %迭代上限1000 if delta<=teminate break; else gk1=gk;fk1=fk;(3)确定步长 %利用Wolfe-Powell搜索计算步长 alphak,fk,gk,wfnum,wgnum=wolfe2(n,m,xk,dk,fk1,gk1,nprob); %利用Wolfe-Powell搜索计算步长 (4)计算 fnum=fnum+wfnum; gnum=gnum+wgnum; xk1=xk;xk=xk1+alphak*dk;(5)确定下降方向 fk=objfcn(n,m,xk,nprob); gk=grdfcn(n,m,xk,nprob); temk1=norm(gk1,2); temk=norm(gk,2); dk1=dk; bk=temk2/temk12; dk=-gk+bk*dk1; end k=k+1; delta=norm(gk,2);End(6)无约束问题运算结束后记录所花费时间time=toc;%终止计时if time<=0. t(i,s)=0.0001; else t(i,s)=time;End(8)输出无约束问题的运行结果fprintf('nt%sttt%2dttt%5dtttt%5dttt%5dttt%4fn',filename,n,k,fnum,gnum,time);End非精确Wolf-powell线性搜索function alphak1,fk2,gk2,wfnum,wgnum=wolfe2(n,m,xk,dk,fk,gk,nprob) rho1=0.8;rho2=0.6;sigma1=0.01;sigma2=0.6;% fk1=fk;gk1=gk; wfnum=0;wgnum=0; %step 0 alphak1=1; fk2=objfcn(n,m,xk+alphak1*dk,nprob);wfnum=wfnum+1; gk2=grdfcn(n,m,xk+alphak1*dk,nprob);wgnum=wgnum+1; if fk2-fk1<=sigma1*alphak1*gk1'*dk if gk2'*dk>=sigma2*gk1'*dk return; end %step 0 else %step 1% i=-8; while 1 if i=0 alphak1=rho1i; fk2=objfcn(n,m,xk+alphak1*dk,nprob);wfnum=wfnum+1; end if fk2-fk1<=sigma1*alphak1*gk1'*dk i=i-1; fk2=objfcn(n,m,xk+rho1i*dk,nprob); if fk2-fk1>sigma1*rho1i*gk1'*dk %alphak=alphak1 break; end else i=i+1; end end %alphak1=rho1i %step 1% j=0; while 1 alphak1=rho2j*alphak1; if alphak1=0 break; end gk2=grdfcn(n,m,xk+alphak1*dk,nprob);wgnum=wgnum+1; if gk2'*dk>=sigma2*gk1'*dk %alphak=alphak1; break; end j=j+1; end End(1) 参数设置：Wolf-powell搜索中的；拟牛顿法及共轭梯度法中相关参数:teminate=1.0e-6;factor=0.1numer=1,2,3,4,5,6,7,8,9,11,12,13,14,15,16,18,19,20,21, 22,23,24,25,26,28,29,30,31,32,33,34;(2) 终止准则：Wolf-powell算法终止：当搜索到步长满足（3.1）的第二个不等式时终止程序；拟牛顿法及共轭梯度法算法终止：当时，此处，迭代次数，若迭代次数达到1000，仍无法满足的条件，则退出算法。三、 实验结果及分析运行Matlab程序（见附录）得如下结果：*拟牛顿法results*ProblemDim.Iter.fnumgnumtime*rose 2 327 362 3301.froth 2 202 228 2040.badscp 2 1000 1081 10020.badscb 2 156 219 1590.beale 2 394 395 3950.jensam 2 81 109 840.helix 3 131 212 1360.bard 3 1000 1052 10032.gauss 3 771 772 7721.gulf 3 1000 1001 10011.box 3 262 263 2630.sing 4 1000 1028 10030.wood 4 245 365 2540.kowosb 4 1000 1001 10011.bd 4 100011754601.bigss 6 1000 1014 10026.osb211 1000 1050 10044.watson100 1000 1172 100932.rosex100 427 1651 4954.singx20 1000 1028 10032.pen110 1000 1001 10013.pen210 1000 1001 10012.vardim10 84 148 860.trig 10 62 63 630.bv 10 1000 1001 10011.IE 100 60 61 613.trid 10 31 180 460.band10 28 241 460.lin 10 87 88 880.lin1 10 4 56 60.lin0 10 4 52 60.*FR共轭梯度results*ProblemDim.Iter.fnumgnumtime*rose 2 219 7165 4391.froth 2 234 9462 4691.badscp 2 100096436 200115.badscb 2 1000 200715.beale 2 255 6252 5030.jensam 2 97 3736 1950.helix 3 251 8532 5031.bard 3 362 7269 7193.gauss 3 221 4034 4120.gulf 3 1 13 30.box 3 100041907 19947.sing 4 100030599 20013.wood 4 100056767 20015.kowosb 4 100039371 19876.bd 4 1000 8409317.bigss 6 100028339 19984.osb211 100039071 200116.watson100 100063667 2001351.rosex100 43915429 8795.singx20 100029525 20017.pen110 46011699 5368.pen210 63415400 12553.vardim10 138 4896 2770.trig 10 34 179 510.bv 10 90521683 18114.IE 100 107 1124 20312.trid 10 100028795 20013.band10 1000717831237013.lin 10 44 243 670.lin1 10 84 5070 1690.lin0 10 1000 10278321.运行时间分布图对比：结果分析：总体上可明显看出拟牛顿法比FR共轭梯度法效率要高。虽然对某一特定的无约束问题，可能会出现用拟牛顿法反而没有FR共轭梯度法节约时间，计算量小，迭代次数少的情况。从得出的数据上看，由于34个不同的问题，除个别问题外，拟牛顿法运行时间基本保持在很短的水平，故得出的图像较为平坦，波动较小；而共轭梯度法在处理不同问题时所用的时间较不稳定，且有些问题耗费时间明显较多，故图像浮动较大。故可以得出结论，较FR共轭梯度法，拟牛顿法在解决无约束问题上效率更高，稳定性更好。四、 实验心得在平时上课时有很多求解无约束问题的方法只知道一些大概，许多算法上的细节上缺乏深入的理解，通过这次实验在这方面知识上得到了补充，同时在Matlab软件应用上也有了一定的进步。实验的准备期，通过翻阅课本及网上资源的查找，我逐渐理解了最速下降法、牛顿法、拟牛顿法及共轭梯度法在算法上的区别，各个算法的性质，存在的优点和缺陷。了解到这些方法之间并不是独立的，牛顿法可以说是在最速下降法上的升华，拟牛顿法又是在克服牛顿法的缺陷的基础上建立的，共轭梯度法则是介于最速下降法和牛顿法之间的一种方法。解决无约束问题的算法百样，但是我们现所学的都是基于下降算法的方法，无论是最速下降法，拟牛顿法，或是共轭梯度法，都是主要解决如何寻找更优更快的下降方向，及如何更加合理地取得步长，这一共同点使我在理解算法上得到了很大的帮助。Matlab软件实现过程，由于Matlab软件是这学期夏季学期学习了，虽然一些基本语法还记得，但是在编写程序上还是存在困难，经常遇到报错的情况，通过多次的调试才得以运行。不过正是这个算法转化为程序语言的过程，让我较之前能更熟练得运用Matlab软件，也更加注意算法中的一些细节。附录：%拟牛顿法及FR共轭梯度法 program using Wolfe-Powell search %clc;muk=10;teminate=1.0e-6;factor=0.1;numer=1,2,3,4,5,6,7,8,9,11,12,13,14,15,16,18,19,20,21,22,23,24,25,26,28,29,30,31,32,33,34;no=size(numer);nn=no(:,2);s=1;fprintf('nt*拟牛顿法results*n');fprintf('tProblemttDim.ttIter.tttfnumttgnumtttimen');fprintf('t*n');for i=1:nnnprob=numer(i); n,m,xk,filename=initf(nprob);% 读初始数据xk=factor*xk;bk=eye(n);k=0;tic; %计时开始 fk=objfcn(n,m,xk,nprob);fnum=1;gk=grdfcn(n,m,xk,nprob);gnum=1;delta=norm(gk,2);while k<1000 %迭代上限1000 if delta<=teminate break; else dk=-linsolve(bk+muk*eye(n),gk); gk1=gk;fk1=fk;gkdk=gk'*dk; if gk'*dk>=-1.0e-14 %当dk不是充分下降时采用负梯度为搜索方向 dk=-gk; end %利用Wolfe-Powell搜索计算步长 alphak,fk,gk,wfnum,wgnum=wolfe2(n,m,xk,dk,fk1,gk1,nprob); %利用Wolfe-Powell搜索计算步长 fnum=fnum+wfnum; gnum=gnum+wgnum; xk1=xk;xk=xk1+alphak*dk; fk=objfcn(n,m,xk,nprob); gk=grdfcn(n,m,xk,nprob); if norm(gk,2)<=teminate k=k+1; break; end % Bk update sk=xk-xk1;bks2=sk'*bk*sk;yk=gk-gk1; yksk=yk'*sk; if yksk>0 bks1=bk*sk*sk'*bk; yks=yk*yk'/yksk;bk1=bk; bk=bk1-bk1*sk*sk'*bk1/(sk'*bk1*sk)+yk*yk'/(yk'*sk); end end k=k+1;endtime=toc;%终止计时if time<=0. t(i,s)=0.0001; else t(i,s)=time;endfprintf('nt%stt%2dtt%5dttt%5dtt%5dtt%4fn',filename,n,k,fnum,gnum,time);ends=2;fprintf('nt*FR共轭梯度法results*n');fprintf('tProblemttDim.ttIter.tttfnumttgnumtttimen');fprintf('t*n');for i=1:nnnprob=numer(i); n,m,xk,filename=initf(nprob);% 读初始数据xk=factor*xk;bk=eye(n);k=0;tic; %计时开始 fk=objfcn(n,m,xk,nprob);fnum=1;gk=grdfcn(n,m,xk,nprob);gnum=1;delta=norm(gk,2);dk=-gk;while k<1000 %迭代上限1000 if delta<=teminate break; else gk1=gk;fk1=fk; %利用Wolfe-Powell搜索计算步长 alphak,fk,gk,wfnum,wgnum=wolfe2(n,m,xk,dk,fk1,gk1,nprob); %利用Wolfe-Powell搜索计算步长 fnum=fnum+wfnum; gnum=gnum+wgnum; xk1=xk;xk=xk1+alphak*dk; fk=objfcn(n,m,xk,nprob); gk=grdfcn(n,m,xk,nprob); temk1=norm(gk1,2); temk=norm(gk,2); dk1=dk; bk=temk2/temk12; dk=-gk+bk*dk1; end k=k+1; delta=norm(gk,2);endtime=toc;%终止计时if time<=0. t(i,s)=0.0001; else t(i,s)=time;endfprintf('nt%stt%2dtt%5dttt%5dtt%5dtt%4fn',filename,n,k,fnum,gnum,time);end % clc;%作图开始ZT=0.01;Tau=10;Mp=size(t,1);%测试问题数目%Ms=size(t,2);%求解方法数目%r=zeros(Mp,Ms);for i=1:Mp mintp=min(t(i,1:Ms);%在所有求解器中对每一测试问题提取最小CPU时间% if mintp=0 mintp=0.001; end for s=1:Ms r(i,s)=t(i,s)/mintp; endendrho=zeros(Ms);hold on; for s=1:Ms numbers=zeros(1,(Tau)/ZT); for tau=0.00001:ZT:Tau for i=1:Mp if r(i,s)<=tau k=ceil(tau-1)/ZT); numbers(k)=numbers(k)+1; end end end switch s case 1 % 拟牛顿法 tau=1.00001:ZT:Tau; k=ceil(tau-1)/ZT); plot(tau,numbers(k)/Mp,':b'); case 2 % FR共轭梯度法 tau=1.00001:ZT:Tau; k=ceil(tau-1)/ZT);
