1、数学实验报告 实验 6 无约束优化1实验 6 无约束优化分 1 黄浩 2011011743一、 实验目的1. 掌握用 MATLAB 优化工具箱的基本用法，对不同算法进行初步分析、比较。2. 练习用无约束优化方法建立和求解实际问题模型（包括非线性最小二乘拟合）。二、 实验内容1. 数学实验第二版（问题 2.1）问题叙述：取不同的初值计算非线性规划： min=(12)2(11)2112(11)52尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法（搜索方向、步长搜索、数值梯度与分析梯度等）的结果进行分析、比较。实验过程：首先绘制这个函数的三维图形以及等高线（程序见四.1） ，结果如下：数学
2、实验报告 实验 6 无约束优化2通过观察这两幅图，可以得到，x2 确定时，x1 越负，函数值越大，x1 确定时，x2 绝对值越大，函数值越大。但对于 x1 正向偏离 0 的情况，并没有很好的反映，于是扩大绘图范围，做出下图（程序见四.2）：数学实验报告 实验 6 无约束优化3由上面两幅图可见，方程像是一个四角被捏起的花布，而且 z 的最小值为0。因此只要求解该方程的零点，即得到了方程的局部极小点，且若将原方程变形为：=(12)2(11)412(11)42我们容易发现，该方程的零点为：x2=0 或 x1=0 或 x1=1 或 2(11)4=1在求解零点之前，先针对一个零点，不妨用 x1=1，x2
3、=1，分析不同算法的优劣。在 matlab 的无约束优化中，可以使用 fminumc 和 fminsearch 两种函数，搜索方向的算法有 BFGS 公式、DFP 公式和最速下降法三种（书中还提到的Gill-Murray 公式在 matlab 中已经不再使用） ，步长的一维搜索有混合二次三次多项式插值和三次多项式插值两种方法，另外，在求解函数梯度是也有数值方法和分析方法两种。在对上述四类算法因素进行分析时，我们采用控制变量法，每次只保持一种或两种算法因素改变，分析它的精度及效率。（一） 分析 fminumc 与 fminsearch 两种方法的精度及效率选择初值为 x1=0.8,x2=0.8，
4、使用 fminunc 和 fminsearch 的默认算法及控制参数，输出结果如下（程序见四.3、四.4）：数学实验报告 实验 6 无约束优化4最优解 x1 最优解 x2 最优值 z 迭代次数 z 调用次数fminunc 9.9400e-001 7.9238e-001 8.0435e-010 12 39fminsearch 1.0000e+000 5.8842e-001 9.5246e-057 67 123因为精确解为 x1=1，z=0，我们便可以比较出不同算法的精度。由上表可见，fminunc 的迭代次数和目标函数的调用次数较少，效率较高，但最优值的精度较低，与之相反的是，fminsearc
5、h 的迭代次数和目标函数的调用次数较多，效率比较低，然而精度很高。（二） 比较不同搜索方向和不同步长的算法精度及效率使用 fminunc 函数，选择初值为 x1=0.6,x2=0.8，分别用 BFGS 公式、DFP公式和最速下降法三种搜索方向，混合二次三次多项式插值和三次多项式插值两种步长算法，共 2*3=6 种方式进行运算，结果如下（程序见四.5）：搜索方向 步长算法 最优解 x1 最优解 x2 最优值 z 迭代次数 z 调用次数BFGS 1.0072e+000 7.6715e-001 1.5634e-009 7 27DFP 1.0074e+000 7.5404e-001 1.7336e-0
6、09 5 21最速下降法混合二次三次多项式插值 9.9290e-001 7.2793e-001 1.3266e-009 39 255BFGS 1.0072e+000 7.6715e-001 1.5634e-009 7 27DFP 1.0074e+000 7.5404e-001 1.7336e-009 5 21最速下降法三次多项式插值9.9290e-001 7.2793e-001 1.3266e-009 39 255显然，精确解应为 x1=1，z=0.由上表可知，计算效率：DFPBFGS最速下降法，而且 DFP 和 BFGS 的迭代次数近似，最速下降法的迭代次数远远大于前两者。这是因为最速下降法
7、的收敛阶数为 1，开始时收敛非常快，但当接近最优解时，收敛就变慢了。而前两种算法均利用到了二阶导数项，收敛阶数较高，因此总体上收敛速率较快，计算效率更高。对于解的精度，三者近似相同，没有很大的区别。 此外，通过对比可以看出，对于混合二次三次多项式插值和三次多项式插值两种步长算法，他们的最优解以及迭代次数都是一样的。但事实上不是二者没有区别，我通过翻阅 matlab 的 warning 消息，得知 LineSearchType 的控制参数已经不存在了，因而 matlab 只使用了某一种步长算法。（三） 比较数值方法和分析方法的精度及效率首先编写含梯度的函数 M 文件（程序见四.6） ，然后选择初
8、值数学实验报告 实验 6 无约束优化5x1=0.6，x2=0.8，分别使用数值方法和分析方法进行计算（程序见四.7） ，结果如下：搜索方向 算法 最优解 x1 最优解 x2 最优值 迭代次数 z 调用次数BFGS 1.00720E+00 7.67150E-01 1.56340E-09 7 27DFP 1.00740E+00 7.54040E-01 1.73360E-09 5 21最速下降法数值方法9.92900E-01 7.27930E-01 1.32660E-09 39 255BFGS 1.00720E+00 7.67150E-01 1.56340E-09 7 9DFP 1.00740E+0
9、0 7.54040E-01 1.73360E-09 5 7最速下降法分析方法9.92900E-01 7.27930E-01 1.32660E-09 39 85由上表可知，分析方法的目标函数调用次数少于数值方法，在进行大规模计算的时候，分析方法不失为数值方法的一种替代。但从本例来看，分析方法对于提高精度的效果不明显，而且对迭代次数也没有减少。在进行完上述算法的比较之后，我们开始逐个检验之前得到的极小值点，选定 fminunc 函数和 BFGS 搜索方向进行分析, 取不同的初值（程序见四.8） ，结果如下：初值 x01 初值 x02 期望解 X1 期望解 X2 最优解 x1 最优解 x2 最优值0
10、.3 0.3 0 0 2.2336e-001 2.8437e-006 1.4678e-0130.2 -2 0 - -7.5875e-009 -1.3095e+000 5.2657e-0161.2 4 1 - 1.0026e+000 3.9964e+000 7.2181e-0100.3 18 0.5 16 1.0088e+000 1.7980e+001 1.9765e-0060.45 15.6 0.5 16 4.9682e-001 1.5599e+001 2.6331e-014在第四步实验（即上表第四行）中，作者本来希望初值（0.3,18）能收敛到期望值（0.5,16） ，然而却最终得到的最优解
11、为（1,18） ，虽然也为极值点，但搜索方向却没有按照预期进行。于是进行了第五步的实验，将初值与期望值的偏差进一步缩小，结果获得了期望的最优解。由上表的最优值，我们可以判定，当 x2=0 或 x1=0 或 x1=1 或时，目标函数值为 0，又因为目标函数是非负连续的，因而零点2(11)4=1即为局部极小值，也为全局极小值。但这只是局部极小值的充分条件，仅凭此并不能找到所有的局部极小值。数学实验报告 实验 6 无约束优化6然而，通过对图像的观察，可以直观地判断出零点和极小值点是等价的，至于更分析上的探讨，由于本人数学能力有限，无法做出全面的证明。得出结论：该函数的极小值点为：x2=0 或 x1=
12、0 或 x1=1 或 ，极小值点2(11)4=1即为零点，因而极小值为 0.对于不同算法的优劣，有如下四点结论：a) fminunc 的迭代次数和目标函数的调用次数较少，效率较高，但最优值的精度较低，与之相反的是，fminsearch 的迭代次数和目标函数的调用次数较多，效率比较低，然而精度很高。b) 对于 DFP、BFGS 和最速下降法这三种搜索方向，计算效率：DFPBFGS最速下降法，而三者的计算精度在数量级上是近似的。c) 对于混合二次三次多项式插值和三次多项式插值两种步长算法，由于matlab 中已经没有了步长的设置参数，因此无法进行比较。d) 对于数值方法和分析方法计算梯度，结论是分
13、析方法可以大幅减少目标函数的调用次数，但对于提高精度的贡献很小，迭代次数也没有区别。2. 数学实验第二版（问题 6）问题描述：有一组数据(t i,yi)(i=1,2,33)，其中 ti=10(i-1)，y i由表中给出。现在用这组数据拟合函数： (,)=1+24+35其中参数 x，初值可选为(0.5,1.5,-1,0.01,0.02)，用 GN 和 LM 两种方法求解。对 yi 作一扰动，即 yi+ei，e i 为(-0.05,0.05)内的随机数，观察并分析迭代收敛是否会变慢。 1 0.844 12 0.718 23 0.478 2 0.908 13 0.685 24 0.467 3 0.9
14、32 14 0.658 25 0.457 4 0.936 15 0.628 26 0.448 5 0.925 16 0.603 27 0.438 6 0.908 17 0.580 28 0.431 7 0.881 18 0.558 29 0.424 数学实验报告 实验 6 无约束优化78 0.850 19 0.538 30 0.420 9 0.818 20 0.522 31 0.414 10 0.784 21 0.506 32 0.411 11 0.751 22 0.490 33 0.406 实验过程：在 matlab 中，实现最小二乘拟合有 lsqnonlin 和 lsqcurvefit，鉴
15、于后者的代码更清晰简洁，因而我使用了 lsqcurvefit 来进行拟合。又因为在此版本的 matlab 中，已经不能设置 LevenbergMarquardt 的控制参数，即不能手动设置非线性最小二乘拟合的下降方向，因此下面我只给出了LM 法的拟合结果（程序见四.9、四.10）：x1 x2 x3 x4 x5 迭代数 z 调用次数 norm0.37541 1.9358 -1.4647 0.01287 0.02212 8 54 5.4649e-005因此，对于函数 ，五个参数的取值分别为(,)=1+24+350.37541、1.9358、-1.4647、0.01287、0.02212 时，与上述
16、数据点的差别最小，拟合度最高。误差平方和为 5.4649e-005然后，对 yi 作一扰动，即 yi+ei，e i 为(-0.05,0.05)内的随机数，并重复四次（程序见四.11）所得结果如下：序号 x1 x2 x3 x4 x5 迭代数 z 调用次数 norm1 0.401 2.399 -1.996 0.015 0.023 9 60 2.2390e-0022 0.342 1.307 -0.817 0.010 0.025 8 54 2.2705e-0023 0.407 25.920 -25.500 0.018 0.018 50 306 3.0588e-0024 0.367 1.540 -1.1
17、04 0.012 0.026 11 72 2.1883e-002由上表可以看出，扰动后输出的误差平方和 norm 数量级在 0.01 左右，远远大于未扰动之前的误差。因此，当一个随机小扰动加到 y 上之后，新的数据点一般会产生较大的偏离，致使数据点与函数的吻合性变差了。不仅如此，当加入扰动后，迭代数与目标函数调用次数也有或多或少的增加，这说明收敛速度变慢了，这也是由于数据点拟合度不高的必然结果。得出结论：对于函数 ，五个参数的取值分别为(,)=1+24+35数学实验报告 实验 6 无约束优化80.37541、1.9358、-1.4647、0.01287、0.02212 时，拟合度最高。而且在加
18、入扰动后，拟合度和精度大大下降，收敛速度也有所变慢。三、 实验总结本次实验是利用 matlab 优化工具箱进行无约束优化，并进行多元函数极小值的分析及非线性含参量函数的最小二乘拟合。这一部分无论在基本原理上还是代码上都比之前的实验更有难度，如果说之前的实验是灰箱的话，那么我认为本次实验完全是“半黑”箱，各种搜索方向和步长算法完全是照搬代码，而对其内部的原理不甚清楚，这就导致在得出实验结果之后，只能进行表面上的分析而不能从内部原理来解释各种算法的优劣。当然，对于今后可能的工程应用而言，数学的工具性高于它的理论性，通过这次实验，我也基本掌握了函数极值的基本求解思路及非线性最小二乘的方法，对算法和精
19、度控制的有关参数也能够熟练应用，达到了基本要求。四、 程序清单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(
20、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 的优劣数学实验报告 实验 6 无约束优化9x0=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=
21、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,s
22、teepdesc,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. 第一题设
23、置分析方法的函数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 nargout1g(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)
24、*(1-x(1)5);end7. 第一题比较数值方法和分析方法计算梯度数学实验报告 实验 6 无约束优化10x0=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=optim
25、set(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
26、,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.
27、 第一题选定 fminunc 函数和 BFGS 搜索方向, 取不同的初值以寻找极小值点x0=10,6; %x0可以变化%format shorteopt=optimset(LargeScale,off,MaxFunEvals,1000,HessUpdate,bfgs);x1,z1,ef1,out1=fminunc(profun1,x0,opt)9. 第二题设置非线性最小二乘拟合的函数function f=profun2(x,t)f=x(1)+x(2)*exp(-x(4)*t)+x(3)*exp(-x(5)*t);end10. 第二题无扰动时，最小二乘拟合i=0:32;t=10*i;y=.844
28、 .908 .932 .936 .925 .908 .881 .850 .818 .784 .751 .718 .685 .658 .628 .603 .580 .558 .538 .522 .506 .490 .478 .467 .457 .448 .438 .431 .424 .420 .414 .411 .406;x0=0.5,1.5,-1,0.01,0.02;数学实验报告 实验 6 无约束优化11x,norm,res,ef,out,lam,jac=lsqcurvefit(profun2,x0,t,y)11. 第二题有扰动时，最小二乘拟合i=0:32;t=10*i;e=0.1*rand(1,33)-0.05;y=.844 .908 .932 .936 .925 .908 .881 .850 .818 .784 .751 .718 .685 .658 .628 .603 .580 .558 .538 .522 .506 .490 .478 .467 .457 .448 .438 .431 .424 .420 .414 .411 .406;x0=0.5,1.5,-1,0.01,0.02;y=y+e;x,norm,res,ef,out,lam,jac=lsqcurvefit(profun2,x0,t,y)
