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