數(shù)學(xué)實(shí)驗(yàn)——無約束優(yōu)化.docx

數(shù)學(xué)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)6 無約束優(yōu)化實(shí)驗(yàn)6 無約束優(yōu)化分1 黃浩 2011011743一、 實(shí)驗(yàn)?zāi)康?. 掌握用MATLAB 優(yōu)化工具箱的基本用法，對(duì)不同算法進(jìn)行初步分析、比較。2. 練習(xí)用無約束優(yōu)化方法建立和求解實(shí)際問題模型（包括非線性最小二乘擬合）。二、 實(shí)驗(yàn)內(nèi)容1. 數(shù)學(xué)實(shí)驗(yàn)第二版（問題2.1）問題敘述：取不同的初值計(jì)算非線性規(guī)劃：minz=x1x22(1-x1)21-x1-x2(1-x1)52盡可能求出所有局部極小點(diǎn)，進(jìn)而找出全局極小點(diǎn)，并對(duì)不同算法（搜索方向、步長(zhǎng)搜索、數(shù)值梯度與分析梯度等）的結(jié)果進(jìn)行分析、比較。實(shí)驗(yàn)過程：首先繪制這個(gè)函數(shù)的三維圖形以及等高線（程序見四.1），結(jié)果如下： 通過觀察這兩幅圖，可以得到，x2確定時(shí)，x1越負(fù)，函數(shù)值越大，x1確定時(shí)，x2絕對(duì)值越大，函數(shù)值越大。但對(duì)于x1正向偏離0的情況，并沒有很好的反映，于是擴(kuò)大繪圖范圍，做出下圖（程序見四.2）：由上面兩幅圖可見，方程像是一個(gè)四角被捏起的花布，而且z的最小值為0。因此只要求解該方程的零點(diǎn)，即得到了方程的局部極小點(diǎn)，且若將原方程變形為：z=x1x22(1-x1)41-x2(1-x1)42我們?nèi)菀装l(fā)現(xiàn)，該方程的零點(diǎn)為：x2=0或x1=0或x1=1或x2(1-x1)4=1在求解零點(diǎn)之前，先針對(duì)一個(gè)零點(diǎn)，不妨用x1=1，x2=1，分析不同算法的優(yōu)劣。在matlab的無約束優(yōu)化中，可以使用fminumc和fminsearch兩種函數(shù)，搜索方向的算法有BFGS公式、DFP公式和最速下降法三種（書中還提到的Gill-Murray公式在matlab中已經(jīng)不再使用），步長(zhǎng)的一維搜索有混合二次三次多項(xiàng)式插值和三次多項(xiàng)式插值兩種方法，另外，在求解函數(shù)梯度是也有數(shù)值方法和分析方法兩種。在對(duì)上述四類算法因素進(jìn)行分析時(shí)，我們采用控制變量法，每次只保持一種或兩種算法因素改變，分析它的精度及效率。（一） 分析fminumc與fminsearch兩種方法的精度及效率選擇初值為x1=0.8,x2=0.8，使用fminunc和fminsearch的默認(rèn)算法及控制參數(shù)，輸出結(jié)果如下（程序見四.3、四.4）：最優(yōu)解x1最優(yōu)解x2最優(yōu)值z(mì)迭代次數(shù)z調(diào)用次數(shù)fminunc9.9400e-0017.9238e-0018.0435e-0101239fminsearch1.0000e+0005.8842e-0019.5246e-05767123 因?yàn)榫_解為x1=1，z=0，我們便可以比較出不同算法的精度。由上表可見，fminunc的迭代次數(shù)和目標(biāo)函數(shù)的調(diào)用次數(shù)較少，效率較高，但最優(yōu)值的精度較低，與之相反的是，fminsearch的迭代次數(shù)和目標(biāo)函數(shù)的調(diào)用次數(shù)較多，效率比較低，然而精度很高。（二） 比較不同搜索方向和不同步長(zhǎng)的算法精度及效率使用fminunc函數(shù)，選擇初值為x1=0.6,x2=0.8，分別用BFGS公式、DFP公式和最速下降法三種搜索方向，混合二次三次多項(xiàng)式插值和三次多項(xiàng)式插值兩種步長(zhǎng)算法，共2*3=6種方式進(jìn)行運(yùn)算，結(jié)果如下（程序見四.5）：搜索方向步長(zhǎng)算法最優(yōu)解x1最優(yōu)解x2最優(yōu)值z(mì)迭代次數(shù)z調(diào)用次數(shù)BFGS混合二次三次多項(xiàng)式插值1.0072e+0007.6715e-0011.5634e-009727DFP1.0074e+0007.5404e-0011.7336e-009521最速下降法9.9290e-0017.2793e-0011.3266e-00939255BFGS三次多項(xiàng)式插值1.0072e+0007.6715e-0011.5634e-009727DFP1.0074e+0007.5404e-0011.7336e-009521最速下降法9.9290e-0017.2793e-0011.3266e-00939255顯然，精確解應(yīng)為x1=1，z=0.由上表可知，計(jì)算效率：DFPBFGS最速下降法，而且DFP和BFGS的迭代次數(shù)近似，最速下降法的迭代次數(shù)遠(yuǎn)遠(yuǎn)大于前兩者。這是因?yàn)樽钏傧陆捣ǖ氖諗侩A數(shù)為1，開始時(shí)收斂非?？?，但當(dāng)接近最優(yōu)解時(shí)，收斂就變慢了。而前兩種算法均利用到了二階導(dǎo)數(shù)項(xiàng)，收斂階數(shù)較高，因此總體上收斂速率較快，計(jì)算效率更高。對(duì)于解的精度，三者近似相同，沒有很大的區(qū)別。 此外，通過對(duì)比可以看出，對(duì)于混合二次三次多項(xiàng)式插值和三次多項(xiàng)式插值兩種步長(zhǎng)算法，他們的最優(yōu)解以及迭代次數(shù)都是一樣的。但事實(shí)上不是二者沒有區(qū)別，我通過翻閱matlab的warning消息，得知LineSearchType的控制參數(shù)已經(jīng)不存在了，因而matlab只使用了某一種步長(zhǎng)算法。（三） 比較數(shù)值方法和分析方法的精度及效率 首先編寫含梯度的函數(shù)M文件（程序見四.6），然后選擇初值x1=0.6，x2=0.8，分別使用數(shù)值方法和分析方法進(jìn)行計(jì)算（程序見四.7），結(jié)果如下：搜索方向算法最優(yōu)解x1最優(yōu)解x2最優(yōu)值迭代次數(shù)z調(diào)用次數(shù)BFGS數(shù)值方法1.00720E+007.67150E-011.56340E-09727DFP1.00740E+007.54040E-011.73360E-09521最速下降法9.92900E-017.27930E-011.32660E-0939255BFGS分析方法1.00720E+007.67150E-011.56340E-0979DFP1.00740E+007.54040E-011.73360E-0957最速下降法9.92900E-017.27930E-011.32660E-093985由上表可知，分析方法的目標(biāo)函數(shù)調(diào)用次數(shù)少于數(shù)值方法，在進(jìn)行大規(guī)模計(jì)算的時(shí)候，分析方法不失為數(shù)值方法的一種替代。但從本例來看，分析方法對(duì)于提高精度的效果不明顯，而且對(duì)迭代次數(shù)也沒有減少。在進(jìn)行完上述算法的比較之后，我們開始逐個(gè)檢驗(yàn)之前得到的極小值點(diǎn)，選定fminunc函數(shù)和BFGS搜索方向進(jìn)行分析, 取不同的初值（程序見四.8），結(jié)果如下：初值x01初值x02期望解X1期望解X2最優(yōu)解x1最優(yōu)解x2最優(yōu)值0.30.3002.2336e-0012.8437e-0061.4678e-0130.2-20-7.5875e-009-1.3095e+0005.2657e-0161.241-1.0026e+0003.9964e+0007.2181e-0100.3180.5161.0088e+0001.7980e+0011.9765e-0060.4515.60.5164.9682e-0011.5599e+0012.6331e-014在第四步實(shí)驗(yàn)（即上表第四行）中，作者本來希望初值（0.3,18）能收斂到期望值（0.5,16），然而卻最終得到的最優(yōu)解為（1,18），雖然也為極值點(diǎn)，但搜索方向卻沒有按照預(yù)期進(jìn)行。于是進(jìn)行了第五步的實(shí)驗(yàn)，將初值與期望值的偏差進(jìn)一步縮小，結(jié)果獲得了期望的最優(yōu)解。由上表的最優(yōu)值，我們可以判定，當(dāng)x2=0或x1=0或x1=1或x2(1-x1)4=1時(shí)，目標(biāo)函數(shù)值為0，又因?yàn)槟繕?biāo)函數(shù)是非負(fù)連續(xù)的，因而零點(diǎn)即為局部極小值，也為全局極小值。但這只是局部極小值的充分條件，僅憑此并不能找到所有的局部極小值。然而，通過對(duì)圖像的觀察，可以直觀地判斷出零點(diǎn)和極小值點(diǎn)是等價(jià)的，至于更分析上的探討，由于本人數(shù)學(xué)能力有限，無法做出全面的證明。得出結(jié)論：該函數(shù)的極小值點(diǎn)為：x2=0或x1=0或x1=1或x2(1-x1)4=1，極小值點(diǎn)即為零點(diǎn)，因而極小值為0.對(duì)于不同算法的優(yōu)劣，有如下四點(diǎn)結(jié)論：a) fminunc的迭代次數(shù)和目標(biāo)函數(shù)的調(diào)用次數(shù)較少，效率較高，但最優(yōu)值的精度較低，與之相反的是，fminsearch的迭代次數(shù)和目標(biāo)函數(shù)的調(diào)用次數(shù)較多，效率比較低，然而精度很高。b) 對(duì)于DFP、BFGS和最速下降法這三種搜索方向，計(jì)算效率：DFPBFGS最速下降法，而三者的計(jì)算精度在數(shù)量級(jí)上是近似的。c) 對(duì)于混合二次三次多項(xiàng)式插值和三次多項(xiàng)式插值兩種步長(zhǎng)算法，由于matlab中已經(jīng)沒有了步長(zhǎng)的設(shè)置參數(shù)，因此無法進(jìn)行比較。d) 對(duì)于數(shù)值方法和分析方法計(jì)算梯度，結(jié)論是分析方法可以大幅減少目標(biāo)函數(shù)的調(diào)用次數(shù)，但對(duì)于提高精度的貢獻(xiàn)很小，迭代次數(shù)也沒有區(qū)別。2. 數(shù)學(xué)實(shí)驗(yàn)第二版（問題6）問題描述：有一組數(shù)據(jù)(ti,yi)(i=1,2,33)，其中ti=10(i-1)，yi由表中給出?，F(xiàn)在用這組數(shù)據(jù)擬合函數(shù)：f(x,t)=x1+x2e-x4t+x3e-x5t其中參數(shù)x，初值可選為(0.5,1.5,-1,0.01,0.02)，用GN和LM兩種方法求解。對(duì)yi作一擾動(dòng)，即yi+ei，ei為(-0.05,0.05)內(nèi)的隨機(jī)數(shù)，觀察并分析迭代收斂是否會(huì)變慢。iyiiyiiyi10.844 120.718 230.478 20.908 130.685 240.467 30.932 140.658 250.457 40.936 150.628 260.448 50.925 160.603 270.438 60.908 170.580 280.431 70.881 180.558 290.424 80.850 190.538 300.420 90.818 200.522 310.414 100.784 210.506 320.411 110.751 220.490 330.406 實(shí)驗(yàn)過程：在matlab中，實(shí)現(xiàn)最小二乘擬合有l(wèi)sqnonlin和lsqcurvefit，鑒于后者的代碼更清晰簡(jiǎn)潔，因而我使用了lsqcurvefit來進(jìn)行擬合。又因?yàn)樵诖税姹镜膍atlab中，已經(jīng)不能設(shè)置LevenbergMarquardt的控制參數(shù)，即不能手動(dòng)設(shè)置非線性最小二乘擬合的下降方向，因此下面我只給出了LM法的擬合結(jié)果（程序見四.9、四.10）：x1x2x3x4x5迭代數(shù)z調(diào)用次數(shù)norm0.375411.9358-1.46470.012870.022128545.4649e-005因此，對(duì)于函數(shù)f(x,t)=x1+x2e-x4t+x3e-x5t，五個(gè)參數(shù)的取值分別為0.37541、1.9358、-1.4647、0.01287、0.02212時(shí)，與上述數(shù)據(jù)點(diǎn)的差別最小，擬合度最高。誤差平方和為5.4649e-005然后，對(duì)yi作一擾動(dòng)，即yi+ei，ei為(-0.05,0.05)內(nèi)的隨機(jī)數(shù)，并重復(fù)四次（程序見四.11）所得結(jié)果如下：序號(hào)x1x2x3x4x5迭代數(shù)z調(diào)用次數(shù)norm10.4012.399-1.9960.0150.0239602.2390e-00220.3421.307-0.8170.0100.0258542.2705e-00230.40725.920-25.5000.0180.018503063.0588e-00240.3671.540-1.1040.0120.02611722.1883e-002由上表可以看出，擾動(dòng)后輸出的誤差平方和norm數(shù)量級(jí)在0.01左右，遠(yuǎn)遠(yuǎn)大于未擾動(dòng)之前的誤差。因此，當(dāng)一個(gè)隨機(jī)小擾動(dòng)加到y(tǒng)上之后，新的數(shù)據(jù)點(diǎn)一般會(huì)產(chǎn)生較大的偏離，致使數(shù)據(jù)點(diǎn)與函數(shù)的吻合性變差了。不僅如此，當(dāng)加入擾動(dòng)后，迭代數(shù)與目標(biāo)函數(shù)調(diào)用次數(shù)也有或多或少的增加，這說明收斂速度變慢了，這也是由于數(shù)據(jù)點(diǎn)擬合度不高的必然結(jié)果。得出結(jié)論：對(duì)于函數(shù)f(x,t)=x1+x2e-x4t+x3e-x5t，五個(gè)參數(shù)的取值分別為0.37541、1.9358、-1.4647、0.01287、0.02212時(shí)，擬合度最高。而且在加入擾動(dòng)后，擬合度和精度大大下降，收斂速度也有所變慢。三、 實(shí)驗(yàn)總結(jié)本次實(shí)驗(yàn)是利用matlab優(yōu)化工具箱進(jìn)行無約束優(yōu)化，并進(jìn)行多元函數(shù)極小值的分析及非線性含參量函數(shù)的最小二乘擬合。這一部分無論在基本原理上還是代碼上都比之前的實(shí)驗(yàn)更有難度，如果說之前的實(shí)驗(yàn)是灰箱的話，那么我認(rèn)為本次實(shí)驗(yàn)完全是“半黑”箱，各種搜索方向和步長(zhǎng)算法完全是照搬代碼，而對(duì)其內(nèi)部的原理不甚清楚，這就導(dǎo)致在得出實(shí)驗(yàn)結(jié)果之后，只能進(jìn)行表面上的分析而不能從內(nèi)部原理來解釋各種算法的優(yōu)劣。當(dāng)然，對(duì)于今后可能的工程應(yīng)用而言，數(shù)學(xué)的工具性高于它的理論性，通過這次實(shí)驗(yàn)，我也基本掌握了函數(shù)極值的基本求解思路及非線性最小二乘的方法，對(duì)算法和精度控制的有關(guān)參數(shù)也能夠熟練應(yīng)用，達(dá)到了基本要求。四、 程序清單1. 第一題繪制第一幅三維圖和等高線x1,x2=meshgrid(-2:0.001:1,-1:0.001:1);z=(x1.*x2).2.*(1-x1).2.*(1-x1-x2.*(1-x1).5).2;mesh(x1,x2,z)pause;contour(x1,x2,z,120)2. 第一題繪制第二幅三維圖和等高線x1,x2=meshgrid(-8:0.01:10,-5:0.01:5);z=(x1.*x2).2.*(1-x1).2.*(1-x1-x2.*(1-x1).5).2;mesh(x1,x2,z)pause;contour(x1,x2,z,220)3. 第一題設(shè)置無約束優(yōu)化的函數(shù)function z=profun1(x)z=(x(1)*x(2)2*(1-x(1)2*(1-x(1)-x(2)*(1-x(1)5)2;end4. 第一題比較fminunc和fminsearch的優(yōu)劣x0=0.8,0.8;format shortex1,z1,ef1,out1=fminunc(profun1,x0)x2,z2,ef2,out2=fminsearch(profun1,x0)5. 第一題比較不同搜索方向和不同步長(zhǎng)的算法精度及效率x0=0.6,0.8;format shorteopt=optimset(LargeScale,off,MaxFunEvals,1000);opt1=optimset(opt,HessUpdate,bfgs);opt2=optimset(opt,HessUpdate,dfp);opt3=optimset(opt,HessUpdate,steepdesc);opt4=optimset(opt,HessUpdate,bfgs,LineSearchType,cubicpoly);opt5=optimset(opt,HessUpdate,dfp,LineSearchType,cubicpoly);opt6=optimset(opt,HessUpdate,steepdesc,LineSearchType,cubicpoly);x,fv,exit,out=fminunc(profun1,x0,opt1)pausex,fv,exit,out=fminunc(profun1,x0,opt2)pausex,fv,exit,out=fminunc(profun1,x0,opt3)pausex,fv,exit,out=fminunc(profun1,x0,opt4)pausex,fv,exit,out=fminunc(profun1,x0,opt5)pausex,fv,exit,out=fminunc(profun1,x0,opt6)pause6. 第一題設(shè)置分析方法的函數(shù)function z,g=profun11(x);z=(x(1)*x(2)2*(1-x(1)2*(1-x(1)-x(2)*(1-x(1)5)2;if nargout1 g(1)=2*x(1)*x(2)2*(1-x(1)-x(2)*(1-x(1)5)*(2*x(1)2-3*x(1)+1)*(1-x(1)-x(2)*(1-x(1)5)+(x(1)3-2*x(1)2+x(1)*(5*x(2)*(1-x(1)4-1); g(2)=2*x(1)2*(1-x(1)2*x(2)*(1-x(1)-x(2)*(1-x(1)5)2-2*x(1)2*(1-x(1)7*x(2)2*(1-x(1)-x(2)*(1-x(1)5);end7. 第一題比較數(shù)值方法和分析方法計(jì)算梯度x0=0.6,0.8;format shorteopt=optimset(LargeScale,off,MaxFunEvals,1000);opt1=optimset(opt,HessUpdate,bfgs);opt2=optimset(opt,HessUpdate,dfp);opt3=optimset(opt,HessUpdate,steepdesc);opt=optimset(LargeScale,off,MaxFunEvals,1000,GradObj,on);opt4=optimset(opt,HessUpdate,bfgs);opt5=optimset(opt,HessUpdate,dfp);opt6=optimset(opt,HessUpdate,steepdesc);x1,fv1,exit,out1=fminunc(profun1,x0,opt1)x2,fv2,exit,out2=fminunc(profun1,x0,opt2)x3,fv3,exit,out3=fminunc(profun1,x0,opt3)x4,fv4,exit,out4=fminunc(profun11,x0,opt4)x5,fv5,exit,out5=fminunc(profun11,x0,opt5)x6,fv6,exit,out6=fminunc(profun11,x0,opt6)x=x1;x2;x3;x4;x5;x6fv=fv1;fv2;fv3;fv4;fv5;fv6iterations=out1.iterations;out2.iterations;out3.iterations;out4.iterations;out5.iterations;out6.iterationsfunccount=out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount8. 第一題選定fminunc函數(shù)和BFGS搜索