數學實驗報告實驗6 無約束優化實驗6 無約束優化分1 黃浩 2011011743一、 實驗目的1. 掌握用MATLAB 優化工具箱的基本用法，對不同算法進行初步分析、比較。2. 練習用無約束優化方法建立和求解實際問題模型（包括非線性最小二乘擬合）。二、 實驗內容1. 數學實驗第二版（問題2.1）問題敘述：取不同的初值計算非線性規劃：minz=x1x22(1-x1)21-x1-x2(1-x1)52盡可能求出所有局部極小點，進而找出全局極小點，并對不同算法（搜索方向、步長搜索、數值梯度與分析梯度等）的結果進行分析、比較。實驗過程：首先繪制這個函數的三維圖形以及等高線（程序見四.1），結果如下： 通過觀察這兩幅圖，可以得到，x2確定時，x1越負，函數值越大，x1確定時，x2絕對值越大，函數值越大。但對于x1正向偏離0的情況，并沒有很好的反映，于是擴大繪圖范圍，做出下圖（程序見四.2）：由上面兩幅圖可見，方程像是一個四角被捏起的花布，而且z的最小值為0。因此只要求解該方程的零點，即得到了方程的局部極小點，且若將原方程變形為：z=x1x22(1-x1)41-x2(1-x1)42我們容易發現，該方程的零點為：x2=0或x1=0或x1=1或x2(1-x1)4=1在求解零點之前，先針對一個零點，不妨用x1=1，x2=1，分析不同算法的優劣。在matlab的無約束優化中，可以使用fminumc和fminsearch兩種函數，搜索方向的算法有BFGS公式、DFP公式和最速下降法三種（書中還提到的Gill-Murray公式在matlab中已經不再使用），步長的一維搜索有混合二次三次多項式插值和三次多項式插值兩種方法，另外，在求解函數梯度是也有數值方法和分析方法兩種。在對上述四類算法因素進行分析時，我們采用控制變量法，每次只保持一種或兩種算法因素改變，分析它的精度及效率。（一） 分析fminumc與fminsearch兩種方法的精度及效率選擇初值為x1=0.8,x2=0.8，使用fminunc和fminsearch的默認算法及控制參數，輸出結果如下（程序見四.3、四.4）：最優解x1最優解x2最優值z迭代次數z調用次數fminunc9.9400e-0017.9238e-0018.0435e-0101239fminsearch1.0000e+0005.8842e-0019.5246e-05767123 因為精確解為x1=1，z=0，我們便可以比較出不同算法的精度。由上表可見，fminunc的迭代次數和目標函數的調用次數較少，效率較高，但最優值的精度較低，與之相反的是，fminsearch的迭代次數和目標函數的調用次數較多，效率比較低，然而精度很高。（二） 比較不同搜索方向和不同步長的算法精度及效率使用fminunc函數，選擇初值為x1=0.6,x2=0.8，分別用BFGS公式、DFP公式和最速下降法三種搜索方向，混合二次三次多項式插值和三次多項式插值兩種步長算法，共2*3=6種方式進行運算，結果如下（程序見四.5）：搜索方向步長算法最優解x1最優解x2最優值z迭代次數z調用次數BFGS混合二次三次多項式插值1.0072e+0007.6715e-0011.5634e-009727DFP1.0074e+0007.5404e-0011.7336e-009521最速下降法9.9290e-0017.2793e-0011.3266e-00939255BFGS三次多項式插值1.0072e+0007.6715e-0011.5634e-009727DFP1.0074e+0007.5404e-0011.7336e-009521最速下降法9.9290e-0017.2793e-0011.3266e-00939255顯然，精確解應為x1=1，z=0.由上表可知，計算效率：DFPBFGS最速下降法，而且DFP和BFGS的迭代次數近似，最速下降法的迭代次數遠遠大于前兩者。這是因為最速下降法的收斂階數為1，開始時收斂非常快，但當接近最優解時，收斂就變慢了。而前兩種算法均利用到了二階導數項，收斂階數較高，因此總體上收斂速率較快，計算效率更高。對于解的精度，三者近似相同，沒有很大的區別。 此外，通過對比可以看出，對于混合二次三次多項式插值和三次多項式插值兩種步長算法，他們的最優解以及迭代次數都是一樣的。但事實上不是二者沒有區別，我通過翻閱matlab的warning消息，得知LineSearchType的控制參數已經不存在了，因而matlab只使用了某一種步長算法。（三） 比較數值方法和分析方法的精度及效率 首先編寫含梯度的函數M文件（程序見四.6），然后選擇初值x1=0.6，x2=0.8，分別使用數值方法和分析方法進行計算（程序見四.7），結果如下：搜索方向算法最優解x1最優解x2最優值迭代次數z調用次數BFGS數值方法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由上表可知，分析方法的目標函數調用次數少于數值方法，在進行大規模計算的時候，分析方法不失為數值方法的一種替代。但從本例來看，分析方法對于提高精度的效果不明顯，而且對迭代次數也沒有減少。在進行完上述算法的比較之后，我們開始逐個檢驗之前得到的極小值點，選定fminunc函數和BFGS搜索方向進行分析, 取不同的初值（程序見四.8），結果如下：初值x01初值x02期望解X1期望解X2最優解x1最優解x2最優值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在第四步實驗（即上表第四行）中，作者本來希望初值（0.3,18）能收斂到期望值（0.5,16），然而卻最終得到的最優解為（1,18），雖然也為極值點，但搜索方向卻沒有按照預期進行。于是進行了第五步的實驗，將初值與期望值的偏差進一步縮小，結果獲得了期望的最優解。由上表的最優值，我們可以判定，當x2=0或x1=0或x1=1或x2(1-x1)4=1時，目標函數值為0，又因為目標函數是非負連續的，因而零點即為局部極小值，也為全局極小值。但這只是局部極小值的充分條件，僅憑此并不能找到所有的局部極小值。然而，通過對圖像的觀察，可以直觀地判斷出零點和極小值點是等價的，至于更分析上的探討，由于本人數學能力有限，無法做出全面的證明。得出結論：該函數的極小值點為：x2=0或x1=0或x1=1或x2(1-x1)4=1，極小值點即為零點，因而極小值為0.對于不同算法的優劣，有如下四點結論：a) fminunc的迭代次數和目標函數的調用次數較少，效率較高，但最優值的精度較低，與之相反的是，fminsearch的迭代次數和目標函數的調用次數較多，效率比較低，然而精度很高。b) 對于DFP、BFGS和最速下降法這三種搜索方向，計算效率：DFPBFGS最速下降法，而三者的計算精度在數量級上是近似的。c) 對于混合二次三次多項式插值和三次多項式插值兩種步長算法，由于matlab中已經沒有了步長的設置參數，因此無法進行比較。d) 對于數值方法和分析方法計算梯度，結論是分析方法可以大幅減少目標函數的調用次數，但對于提高精度的貢獻很小，迭代次數也沒有區別。2. 數學實驗第二版（問題6）問題描述：有一組數據(ti,yi)(i=1,2,33)，其中ti=10(i-1)，yi由表中給出。現在用這組數據擬合函數：f(x,t)=x1+x2e-x4t+x3e-x5t其中參數x，初值可選為(0.5,1.5,-1,0.01,0.02)，用GN和LM兩種方法求解。對yi作一擾動，即yi+ei，ei為(-0.05,0.05)內的隨機數，觀察并分析迭代收斂是否會變慢。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 實驗過程：在matlab中，實現最小二乘擬合有lsqnonlin和lsqcurvefit，鑒于后者的代碼更清晰簡潔，因而我使用了lsqcurvefit來進行擬合。又因為在此版本的matlab中，已經不能設置LevenbergMarquardt的控制參數，即不能手動設置非線性最小二乘擬合的下降方向，因此下面我只給出了LM法的擬合結果（程序見四.9、四.10）：x1x2x3x4x5迭代數z調用次數norm0.375411.9358-1.46470.012870.022128545.4649e-005因此，對于函數f(x,t)=x1+x2e-x4t+x3e-x5t，五個參數的取值分別為0.37541、1.9358、-1.4647、0.01287、0.02212時，與上述數據點的差別最小，擬合度最高。誤差平方和為5.4649e-005然后，對yi作一擾動，即yi+ei，ei為(-0.05,0.05)內的隨機數，并重復四次（程序見四.11）所得結果如下：序號x1x2x3x4x5迭代數z調用次數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由上表可以看出，擾動后輸出的誤差平方和norm數量級在0.01左右，遠遠大于未擾動之前的誤差。因此，當一個隨機小擾動加到y上之后，新的數據點一般會產生較大的偏離，致使數據點與函數的吻合性變差了。不僅如此，當加入擾動后，迭代數與目標函數調用次數也有或多或少的增加，這說明收斂速度變慢了，這也是由于數據點擬合度不高的必然結果。得出結論：對于函數f(x,t)=x1+x2e-x4t+x3e-x5t，五個參數的取值分別為0.37541、1.9358、-1.4647、0.01287、0.02212時，擬合度最高。而且在加入擾動后，擬合度和精度大大下降，收斂速度也有所變慢。三、 實驗總結本次實驗是利用matlab優化工具箱進行無約束優化，并進行多元函數極小值的分析及非線性含參量函數的最小二乘擬合。這一部分無論在基本原理上還是代碼上都比之前的實驗更有難度，如果說之前的實驗是灰箱的話，那么我認為本次實驗完全是“半黑”箱，各種搜索方向和步長算法完全是照搬代碼，而對其內部的原理不甚清楚，這就導致在得出實驗結果之后，只能進行表面上的分析而不能從內部原理來解釋各種算法的優劣。當然，對于今后可能的工程應用而言，數學的工具性高于它的理論性，通過這次實驗，我也基本掌握了函數極值的基本求解思路及非線性最小二乘的方法，對算法和精度控制的有關參數也能夠熟練應用，達到了基本要求。四、 程序清單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. 第一題設置無約束優化的函數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的優劣x0=0.8,0.8;format shortex1,z1,ef1,out1=fminunc(profun1,x0)x2,z2,ef2,out2=fminsearch(profun1,x0)5. 第一題比較不同搜索方向和不同步長的算法精度及效率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. 第一題設置分析方法的函數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. 第一題比較數值方法和分析方法計算梯度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函數和BFGS搜索