1、无约束优化一 实验目的1 掌握 MATLAB 优化工具箱的基本用法，对不同算法进行初步分析、比较。2 练习用无约束优化方法建立和求解实际问题的模型( 包括最小二乘拟合 )。二 实验内容1 取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法(搜索方向、搜索步长、数值梯度与分析梯度等) 的结果进行分析、比较。=(12)2(11)2112(11)52理论计算部分：根据多元微积分中求解函数的(局部) 极值点的方法，首先要计算原函数的 梯度以及。对于原函数Hessian矩阵=(12)2(11)2112(11)52， 2首先计算梯度=12，1=22(12)(11)2
2、112(11)522(11)(12)2112(11)52+2(52(11)41)(12)2(11)2，2=21(12)(11)2112(11)522(11)5(12)2(11)2，然后计算 Hessian矩阵2=212212221 222，此处由于计算所得式子非常复杂，所以不再继续计算。理论上可以通过求解非线性方程组确定驻点，然后将每个驻点分别代入 中判断矩阵是否正定，即可以确定该驻()=0 2点是否为局部最优解。但是，通过上面的计算可以得知，要求解的非线性方程组较为复杂，及时解除解，要再带入 判断正定也比较困难，所以此处不进行解析求解。Hessian矩阵观察原式，可以简单的看出，原式一定大于
3、 0，所以 ， 是全局极小点，1=0或 1 2=0通过将最后一个式子因式分解112(11)52=(11)212(11)42，还可以看出， 也是一个局部极小点。(2， 1)用 MATLAB 解决问题：对于目标函数是 2 维以下的情况，不妨先绘制出可能出现局部最优解的区域的函数图象，这样有利于直观的判断函数在这个区域的大致情况。在 MATLAB 中输入以下命令：画出原图像如下：画出等高线如下：从图像上可以看出，由于所选择的范围太大，导致图像中右边大面积的点的值因为太小，所以相对于较大的值来讲全部被认为是 ，使得无法确定零点。0缩小选取的范围，输入命令如下：画出图像如下：画出等高线如下：对于点 ，输
4、入命令如下：(2， 1)画出图像如下：画出等高线如下：从两幅图中可以大致看出，当 这一片区域中，函数值均为 ，而点1=0或 1或 2=0 0处，函数值也为 。(2， 1) 0下面用不同的初值进行计算首先编制函数的 M 文件：在命令栏中输入以下命令：输出结果如下：然后改变初值进行计算，输出结果如下：初值 最优解 函数值 迭代次数 函数调用次数-1， -1 (0.0000， 0.5000) 0 1 6-2， -2(0.0003， 1.2880) 6.2061007 22 720.1， 0.2 (0.0000， 0.1719) 1.1184018 6 243， 2 (1.8782， 1.6809)
5、8.2624011 9 392， 0.9 (2.0248， 0.9066) 1.1459014 5 30从上表中可以得到以下结论，极小值为 ，但是极小点不唯一。当初值选择的不同时，最优0解会向 收敛。但是初值不同，收敛的速度不同，初值越接1=0或 1或 2=0或者点 (2， 1)近最优解，收敛越快。下面用不同的搜索方向和步长搜索进行计算在命令栏中输入以下命令( 用三种搜索方向( ， 和最速下降法 )以及两种步长搜索(混BFGSDFP合二、三次插值和三次插值)：分析梯度只需将命令中的1=(, ,1000)改为1=(, ,1000,)即可。同时函数文件需要作如下改动(此处式子过长不适宜截图 )：然
6、后在命令栏中输入以下命令：输出结果如下表所示：数值梯度搜索方向 步长搜索 最优解 x 最优解 y 函数值 迭代次数 调用次数 3.79570071.7187001 2.9188015 5 21混合二、三次插值 1.90220071.7109001 7.2776016 5 21function y,g=work1(x)y=(x(1)*x(2)2*(1-x(1)2*(1-x(1)-x(2)*(1-x(1)5)2;if nargout1g(1)=2*x(2)*(x(1)*x(2)*(1-x(1)2*(1-x(1)-x(2)*(1-x(1)5)2-2*(1-x(1)*(x(1)*x(2)2*(1-x(
7、1)-x(2)*(1-x(1)5)2+2*(5*x(2)*(1-x(1)4-1)*(x(1)*x(2)2*(1-x(1)2;g(2)=2*x(1)*(x(1)*x(2)*(1-x(1)2*(1-x(1)-x(2)*(1-x(1)5)2-2*(1-x(1)5*(x(1)*x(2)2*(1-x(1)2;end最速下降 1.6709005 1.6971001 5.5436012 18 111 3.79570071.7187001 2.9188015 5 21 1.90220071.7109001 7.2776016 5 21最速下降三次插值1.6709005 1.6971001 5.5436012
8、18 111分析梯度搜索方向 步长搜索 最优解 x 最优解 y 函数值 迭代次数 调用次数 1.47760051.7567001 4.5788012 5 7 5.74710061.7434001 6.8437013 5 7最速下降混合二、三次插值 1.6238005 1.7418001 5.4550012 22 45 1.47760051.7567001 4.5788012 5 7 5.74710061.7434001 6.8437013 5 7最速下降三次插值1.6238005 1.7418001 5.4550012 22 45分析上表，可得到以下结果：1 函数的极小值为 ，分别在 处取得。
9、0 1=0或 1或 2=0或者点 (2， 1)2 当初值不同时，对于得到的最优解、迭代的次数以及函数的调用次数都会有影响。而且当初值不同时，得到的解不同，且得到的都是局部最优解。3 从表中可以看出，使用拟牛顿法的 公式或者 公式时，需要的迭代次数相差不多， 但是使用最速下降法是需要的迭代系数相对就会较多，函数调用次数也会相应增多。4 分析梯度与数值梯度相比，迭代次数相差不大，相差较多的是函数调用次数。而且在本题中，从数值来看，数值梯度相较于分析梯度要较好一些。下面用自己实现的最速下降法和牛顿法来求解最速下降法在命令栏中输入以下内容：其中前面的 是手动求出的梯度， 。p与 q x1(1)=0.5
10、， 2(1)=0.5为 初 值 。输出迭代的结果， ，迭代次数为 次。x1=0.9614， 2=0.2428 161662 有一组数据 其中 ， 由下表给出。现要求用这(， )(=1， 2， 333)， =10(1)组数据拟合函数(， )=1+24+35中的参数 ，初值可选为 ，用 和 两种方法求解。对 (0.5， 1.5， 1， 0， 0.01， 0.02)GNLM作一扰动，即 ， 为 内的随机数，观察并分析迭代收敛是否会变慢。 + (0.05， 0.05) 1 0.844 12 0.718 23 0.4782 0.908 13 0.685 24 0.4673 0.932 14 0.658
11、25 0.4574 0.936 15 0.628 26 0.4485 0.925 16 0.603 27 0.4386 0.908 17 0.580 28 0.4317 0.881 18 0.558 29 0.4248 0.850 19 0.538 30 0.4209 0.818 20 0.522 31 0.41410 0.784 21 0.506 32 0.41111 0.751 22 0.490 33 0.406初步解决：首先编制函数的 M 文件：然后在命令栏中输入以下命令：将所有参数输入到 MATLAB 中。先用 法计算，输入以下命令：LM此时，先不使用大规模算法，输出结果如下：然后再在
12、命令栏中输入以下命令：此时使用大规模算法，输出结果如下：从两个输出结果可以看出，当使用了大规模算法之后，迭代次数和函数调用次数要明显少于没有使用大规模算法的时候。但是二者的解一样。再用 法计算，输入以下命令，同样的，先不使用大规模算法：GN输出结果如下：然后再在命令栏中输入以下命令：此时使用大规模算法，输出结果如下：从输出结果可以看出，使用了大规模算法之后，迭代次数以及函数调用次数确实减少了很多。但是最终算出的结果是一样的。最终算出的结果如下表所示：1 2 3 4 50.37541 1.93585 1.46469 0.01287 0.02212迭代次数的比较：不使用大规模算法 使用大规模算法法
13、LM 35 8法 9 8可以看出， 法的收敛速度不如 法快。LM 下面研究当 发生扰动之后对于迭代收敛速度的影响。(以下全部不使用大规模算法)首先是当 的所有元素均是相等的时候。e，在命令栏中输入以下内容：e=0.4法：LM输出结果如下：法：输出结果如下：，在命令栏中输入以下内容：e=0.1法：LM输出结果如下：法：输出结果如下：，在命令栏中输入以下内容：e=0.1法，输出结果如下：LM法，输出结果如下：，在命令栏中输入以下内容：e=0.4法，输出结果如下：LM法，输出结果如下：通过以上结果可以得到下面的表格：结果的比较e 1 2 3 4 5法LM 0.3354 1.9358 1.46470.
14、0129 0.02210.04法 0.3354 1.9358 1.46470.0129 0.0221法LM 0.3654 1.9358 1.46470.0129 0.02210.01法 0.3654 1.9358 1.46470.0129 0.0221法LM 0.3854 1.9358 1.46470.0129 0.02210.01法 0.3854 1.9358 1.46470.0129 0.02210.04 法LM 0.4154 1.9358 1.46470.0129 0.0221法 0.4154 1.9358 1.46470.0129 0.0221迭代速度的比较e 迭代次数法LM 350.
15、04法 9法LM 350.01法 9法LM 350.01法 9法LM 350.04法 9通过以上的比较可以看出，当 有一定的小扰动时( )，首先发生变y (0.05， 0.05)化的是 ，从表中可以看出，至始至终在变化的只有 ，而其他的解都没有发生变化。其1 1次，迭代次数与函数调用次数没有发生变化，无论扰动有多大，迭代次数和函数调用次数都还是保持原来的值不变，而且依然是 法的迭代次数要大于 法的迭代次数。LM 然后比较 的元素全部是随机数的情况e在命令栏中输入以下内容：LM法输出结果如下：法 在命令栏中输入以下内容：GN输出结果如下：再换另外一组随机数组进行试验输出结果如下：LM法法 输出结
16、果如下：GN再换另外一组随机数组进行试验输出结果如下：LM法法 输出结果如下：GN再换另外一组随机数组进行试验输出结果如下：LM法法 输出结果如下：GN通过以上结果可以得到下面的表格：结果的比较e 1 2 3 4 5法LM 0.3033 1.0629 0.48490.0078 0.02871法 0.3033 1.0629 0.48500.0078 0.0287法LM 0.3418 1.1178 0.65460.0098 0.03422法 0.3419 1.1184 0.65530.0098 0.0341法LM 0.2827 0.9764 0.40900.0072 0.04433法 0.2827 0.9764 0.40900.0072 0.0443法LM 0.3519 5.7454 5.22930.0135 0.01574法 0.3525 19.706319.19080.0142 0.0149迭代速度的比较e 迭代次数法LM 321法 10法LM 82法 14法LM 133法 12法LM 1384法 105当 的值全部都是随机数时，对于结果的变化以及迭代速度的影响都不可预测，结果e可能变大也可能变小，迭代速度可能变快也可能变慢。
