实验7-无约束优化

实验7无约束优化化工系宋劲2011011842第一题【实验描述11取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法(搜索方向、探索步长、数值梯度与分析梯度等)的结果进行分析、比较：minz=(%i%2)2(1—%i)2[l—%i—%2(1—%i)5]2【实验过程11令fiin72函数为/(%i,%2)=(%1%2)2(!-xi)2[l~x1-x2(l-%i)5]2为了直观认识函数，先画出其三维图形，取两自变量均在[-1.5,4]范围内。看z坐标可知此时取的自变量范围太宽，应缩小范围寻找最小值取[x,y]=meshgrid(-0.1:0.01:1.1,-0.1:0.01:1.1)(x为以上xl,y为以上x2)无论从图像上看还是从表达式上看，都很容易得知全局极小点不可能出现在x,y很大的位置。由等高线可以看出最小值在xl=l或0或x2=0的附近。尝试取不同初值找局部极小点，取x0=[・0.5:0.05:1.5,・0.5:0.05:1.5]看这些不同初值对应的局部极小点位置特点对不同算法做对比用分析方法计算梯度再求解【结果分析11直观认识函数，三维视图(x.y]=(-150.1:4-16014)

由图己经看不出最优解的位置，因为范围取太大而使图像太扁，下面缩小了两自变量范围，再作图得：缩小X』范围0.020.0150.010.00501.5-0.5-0.5以下为等高线由等高线可以知道最优解的位置在xl=0、xl=l、x2=O三条线附近。-0.5-0.5从表达式来看，显然有当xl=l或xl=O或x2=0时函数取最小值（因为该函数恒非负故可知f=0必为最小）

求局部极小点，进而求全局极小点取不同的初值看局部极小点分别在哪些位置,这里取X0两个分量均在・0.5至1.5之间，以0.05为间隔取值。此时算法：搜索方向为拟牛顿法的BFGS公式（缺省值），步长探索为混合二、三次插值。做出此时局部极小点的分布图如下：由该图可以较清楚的看出，在x0=[-0.5:0.005:1.5,・0.5:0.005:1.5]时，局部极小点分布在xl=0xl=lx2=0三条直线以及一条曲线上，前三条直线上取到局部极小点己经说明过，也是十分明显的，也同时可以说明这三条线上的点同时还是全局极小点一一因为z是非负的，取到0时必达到极小。主要需要确定这条曲线有什么规律，以及曲线上的点是否是全局极小点。对原程序做小的修改，只取出XK-0.005,且X2X.1的部分点，此时剩下的点均为曲线上的点。对这部分点，取出它们的函数值并做出散点图如下：

由右侧图像可知，曲线上的点对应的函数值均在0附近，非常接近0。由于0一定是全局极小值，所以这条曲线上的点应该是全局极小点，再仔细观察原表达式，可以做出猜想：这条曲线很可能是x.=—^十，因为当此关系被满-QF足时z恒为0。下面将x.=—与该曲线放在一起对比：(If)很明显看出曲线是重合的，可以得出结论：该曲线就是％=序打’且这条曲线上的点均可取到全局极小值。下面尝试取一些较远的初值，取x0=[4:0.05:6,4:0.05:6]看对应的局部极小点与其函数值：6681922122232425262728放大图像看，这些局部极小值中最小的也在200-300之间，肯定都不是全局极小值。故可以得出结论：最优解为xl=O或xl=l或x2=0或x?=——-~o~（If）分析不同算法（搜索方向、探索步长、数值梯度与分析梯度）结果下面用三种搜索方向以及两种步长探索计算，对比结果：不同算法计算结果比较情况搜索方向步长探索最优值X1最优值X2最优值目标函数调用次数1bfgs混合二、二次插值0.5082889.2872e-0081.3027e-016182dip0.5084118.6966e-0081.1417e-016183steepdesc0.5118872.9219e-0051.2698e-0111414bfgs三次插值0.5082889.2872e-0081.3027e-016185dip0.5084118.6966e-0081.1417e-016186steepdesc0.5118872.9219e-0051.2698e-011141由数据表格可知:（1）此时步长探索不影响计算结果，同样的搜索方向在不同步长探索下有相同值（2）最速下降法的效果明显比两个拟牛顿法的效果差，从X2、最优值、调用次数三个指标均可体现：不仅X2与。偏差最大，最优值与0偏差最大，而且调用目标函数次数最多。（3）总的来说，三个方法都比较准确，并且此时都收敛了。下面对比使用两种不同方法得到梯度对结果的影响：以上数据均用数值方法给出近似的梯度，以下数据则用分析方法计算出目标函数梯度：分析方法计算梯度时结果比较情况搜索方向步长探索最优X1最优X2最优值目标函数调用次数1bfgs混合二、二次插值0.5082368.6828e-0081.1389e-01662dip0.5083618.0881e-0089.8771e-01763steepdesc0.5118762.9214e-0051.2695e-011474bfgs三次插值0.5082368.6828e-008L1389e・01665dip0.5083618.0881e-0089.8771e-01766steepdesc0.5118762.9214e-0051.2695e-01147由数据可以看出：此时得到的结果与数值方法得到梯度算出的结果基本相同，误差不大此时调用目标函数次数减少。总的来说分析方法得到梯度会使结果更好一些。【程序清单1】三维画图[x,y]=meshgrid(・L5:0.1:4,・L5:0.1:4);z=((x.*y).A2).*((l-x).A2).*((1-x-y.*((l-x).A5)).A2);mesh(x,y,z)title([x,y]=(-1.5:0.1:4,5:0.1:4))[x,y]=meshgi'id(-0.1:0.01:1.1,-0.1:0.01:1.1);z=((x.*y).A2).*((l-x).A2).*((1-x-y.*((l-x).A5)).A2);mesh(x,y,z)title(缩小x,y范围)pausecontoiir(x,y,z,300)fun72函数functiony=fun72(x)y=((x(l).*x(2)).-2).*((l-x(l)).-2).*((l-x(l)-x(2).*(l-x(l)).A5)).A2;end改变初值求局部极小点n=l;x01=[-0.5:0.05:1.5];x02=[-0.5:0.05:1.5];fori=l:41forj=l:41x0=[x01(i)X)2(j)];optl=optimset(,LargeScale\'off1,MaxFuiiEvals^lOOO);x=fhiiniinc('hin72',x0,optl);kl(n)=x(l);k2(n)=x(2);n=n+l;endendplot(kl,k2,'*)考察曲线11=1;x01=[-0.5:0.05:1.5];xO2=[・O.5:O.O5:L5];fori=l:41forj=l:41x0=[x01(i)X)2(j)];optl=optimset(rLargeScale','off,rMaxFuiiEvals',1000);[x,v]=fininunc(,fiin72#0,optl);ifx(2)>0.1&x(l)<-0.005fiiii(n)=v;kl(n)=x(l);k2(n)=x(2);endn=n+l;endendplot(kl,k2,'+')对比考察xx0=[-l:0.01:0.1];y=(l./(l-xxO)).A4;plot(kl,k2‘b+'xxO,y,T+)算法选择fonnatshortex0=[0.5,0.1];rcasel:bfgsJiybrid23poly'optl=optimset(,LargeScale'/off,*MaxFuiiEvals*,l000/GradObjr,'on');%用分析方法得梯度[xl,vl,exitl,outl]=fhiiniinc('hin72',x0,optl)'…case2:dfj),hybrid2,3poly-—1fopt=optimset(optl,,HessUpdate7dfpr);[x2,v2,exit2,out2]=fininunc(lfiin72',x0,fopt)，・・・case3:gillmiirray,hybrid2,3poly-—rfopt=optimset(optl,,HessUpdate7steepdescr);[x3,v3,exit3,out3]=fiiiinunc(rfiin72',x0,fopt)'…case4:b0s,3rdpoly-—ropt2=optimset(optl/LineSearchiypeVcubicpolyf);[x4,v4,exit4,out4]=fiiiinunc(rfiiii72\x0,opt2)'…case5:dfp,3rdpoly•…'fopt=optimset(opt2,rHessUpdate7dfpr);[x5,v5,exit5,out5]=fiiiinunc(rfiiii72\x0,fopt)'・・・case6:gillmiirray,3rdpoly'fopt=optimset(opt2,rHessUpdate7steepdescr);[x6,v6,exit6,out6]=fiiiinunc(rfiiii72',x0,fopt)Iresultsofsohitions+++'sohitions=[xl;x2;x3;x4;x5;x6];fiinvalues=[vl;v2;v3;v4;v5;v6];iterations=[outl.fimcCount;out2.fiuicCount;out3.fimcCount;out4.fimcCount;out5.fiincCount;out6.fiincCount;];[sohitionsjiinvalues,iterations]计算梯度时函数表达式function[ffg]=fun72(x)f=((x(l).*x(2)).A2).*((l-x(l)).A2).*((l-x(l)-x(2).*(l-x(l)).A5)).A2;ifnargout>lg(l)=2.*x(l).*(x(2).A2).*((l-x(l)).A2).*((l-x(l)-x(2).*((l-x(l)).A5)).A2)-2.*(l-x(l)).*((x(l).*x(2)).A2).*((l-x(l)-x(2).*((l-x(l)).A5)).A2)+2.*((x(l).*x(2)).A2).*((l-x(l)).人2).*(l-x(l)-x(2).*((l-x(l)).A5)).*(5.*x(2).*(x(l).A4)-l)；g(2)=2.*x(2).*(x(l).A2).*((l-x(l)).A2).*((l-x(l)-x(2).*((l-x(l)).A5)).A2)-2.*((x(l).*x(2)).A2).*((l-x(l)).A2).*(1-x(l)-x(2).*((l-x(l)).-5)).*((l-x(l)).-5);endend第二题【实验描述2】有一组数据(t,,yDi=l,2,3……,33.其中ti=10(M),yi由下表给出，现用这些数据拟合函数f(xj)=*+砂F+也eF中的参数X,初值可选为(0.5,1.5,-1,0.01,0.02).用GN和LM两种方法求解。对*作一扰动，即*+%为(・0.05,0.05)内的随机数,观察并分析迭代收敛是否会变慢。.1Yi.1Yi.1Yi10.844120.718230.47820.908130.685240.46730.932140.658250.45740.936150.628260.44850.925160.603270.43860.908170.580280.43170.881180.558290.42480.850190.538300.42090.818200.522310.414100.784210.506320.411110.751220.490330.406【实验过程2]建立函数文fty=fim76(x,t)该题为非线性拟合问题，可用Isqcurvefit命令。题目要求使用LM与GN两种方法求解：参数LevenbergMaiqiiaidt缺省则为LM法，不缺省则为GN法拟合出x后将函数与数据表格中数据在同一图中对比观察拟合情况给y增加(-0.05,0.05)范围内的扰动，此时直接使用LM法求解：求出此时迭代次数、目标函数引用次数，并多次运行程序，观察散点与拟合函数图象的符合程度。【结果分析2】1.LM法求解使用LM求解，输入结果如下：xl=0.37541.9358-1.46470.01290.0221%x取值norml=5.4649e-005%误差平方和resl=0.0026-0.0045-0.00100.00110.00290.00020.00050.0006-0.0006-0.0006-0.0013-0.00100.0007-0.00170.00080.0003-0.00010.00050.0011-0.0006-0.00050.00120.0003-0.0003-0.0006-0.00080.00100.00070.0012-0.00060.0003-0.0013-0.0003%误差向量exitl=3outl=iterations:35%迭代次数fUncCoimt:230%函数调用次数stepsize:2.6301e・005cgiterations:[]firstorderopt:4.0068e-008algoritliin:T^evenberg-Marquardf%LM算法message:[1x460char]由此可得到函数f(x,t)=0.3754+L9358e-0X)129t-14647e_00221t卜面将函数图像与给出的y,t数据点对照考察，看拟合效果如何:由图可见，数据点基本都在函数曲线上。下面再用GN法求解。

2.GN法求解xl=037541.9358-1.46470.01290.0221nonnl=5.4649e-005resl=0.0026-0.0045-0.00100.00110.00290.00020.00050.0006-0.0006■0.0006-0.0013-0.00100.0007-0.00170.00080.0003-0.00010.00050.0011,■0.0006-0.00050.00120.0003-0.0003-0.0006-0.00080.00100.00070.0012-0.00060.0003-0.0013-0.0003exitl=1outl=iterations:9hmcCount:81stepsize:0.9995cgiterations:[]firstorderopt:[]algoritliin:rmedium-scale:Gauss-Newton,line-searchrmessage:[1x147char]可知此时求得的函数f(x,t)与用LM法求得结果是一致的，自然图像也与上图完全一致。但是可注意到，用GN法迭代次数、函数调用次数相比LM法要少得多。所以此时GN法更优，收敛速度更快。3.增加扰动观察规律为了体现随机数的随机性，多次重复试验，并做出图像:第一次第二次

可见此时迭代次数和函数调用次数都随机变化，并不稳定。有时会出现迭代次数比原先更少的情况，也就是说收敛变快；有时乂出现迭代次数远远超过原先的情况，也就是说收敛变慢很多。但增加扰动后拟合的效果肯定会变差，这由图像很容易对比看出来。【程序清单2]构造函数functiony=fun76(x,t)

y=x(l)+x(2)*(exp(-x(4)*t))+x(3)*(exp(-x(5)*t));endLM法x0=[0.5,1.5,-1,0.01,0.02];fori=l:33t(i)=10*(i-l);endc=[0.8440.9080.9320.9360.9250.9080.8810.8500.8180.7840.7510.7180.6850.6580.6280.6030.5800.5580.5380.5220.5060.4900.4780.4670.4570.4480.4380.4310.4240.4200.4140.4110.406];optl=optimset(,LargeScale7off,rMaxFuiiEvals\1000);[xl,norml,resl,exitl,outl]=lsqcurvefit('hm76',x0,t,c,[],[],optl)作图考察x0=[0:0.1:350];y0=0.3754+1.9358*(exp(-0.0129*x0))+(-1.4647)*(exp(-0.0221*x0));plot(t,c,'+'»xO,yO,'')GN法程序基本与LM法一样，只需改一句：optl=optimset(,LargeScale7off,rMaxFuiiEvals\1000,rLevenbergMarquardt\roff);增加扰动x0=[0.5,1.5,-1,0.01,0.02];fori=l:33t(i)=10*(i-l);endc=[0.844+uiiifriid(-0.05,0.05)c=[0.844+uiiifriid(-0.05,0.05)0.908+umfnid(-0.05,0.05)0.936+unifrnd(・0.05,0.05)0.881+iimfind(-0.05,0.05)0.784+iimfind(-0.05,0.05)0.685+unifhid(・0.05,0.05)0.603+unifhid(・0.05,0.05)0.538+unifhid(・0.05,0.05)0.490+unifhid(・0.05,0.05)0.457+unifhid(・0.05,0.05)0.431+iimfind(-0.05,0.05)0.925