数学实验第5次作业_无约束优化.doc

无约束优化一 实验目的1 掌握MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。2 练习用无约束优化方法建立和求解实际问题的模型(包括最小二乘拟合)。二 实验内容1 取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法(搜索方向、搜索步长、数值梯度与分析梯度等)的结果进行分析、比较。minz=x1x22(1-x1)21-x1-x2(1-x1)52理论计算部分：根据多元微积分中求解函数的(局部)极值点的方法，首先要计算原函数的梯度以及Hessian矩阵。对于原函数z=x1x22(1-x1)21-x1-x2(1-x1)52，xR2首先计算梯度z=zx1zx2，zx1=2x2x1x21-x121-x1-x21-x152-21-x1x1x221-x1-x21-x152+25x21-x14-1x1x221-x12，zx2=2x1x1x21-x121-x1-x21-x152-21-x15x1x221-x12，然后计算Hessian矩阵2z=2zx122zx1x22zx2x12zx22，此处由于计算所得式子非常复杂，所以不再继续计算。理论上可以通过求解非线性方程组zx=0确定驻点，然后将每个驻点分别代入2z中判断矩阵是否正定，即可以确定该驻点是否为局部最优解。但是，通过上面的计算可以得知，要求解的非线性方程组较为复杂，及时解除解，要再带入Hessian矩阵判断正定也比较困难，所以此处不进行解析求解。观察原式，可以简单的看出，原式一定大于0，所以x1=0或1，x2=0是全局极小点，通过将最后一个式子因式分解1-x1-x2(1-x1)52=(1-x1)21-x2(1-x1)42，还可以看出，(2，1)也是一个局部极小点。用MATLAB解决问题：对于目标函数是2维以下的情况，不妨先绘制出可能出现局部最优解的区域的函数图象，这样有利于直观的判断函数在这个区域的大致情况。在MATLAB中输入以下命令：画出原图像如下：画出等高线如下：从图像上可以看出，由于所选择的范围太大，导致图像中右边大面积的点的值因为太小，所以相对于较大的值来讲全部被认为是0，使得无法确定零点。缩小选取的范围，输入命令如下：画出图像如下：画出等高线如下：对于点(2，1)，输入命令如下：画出图像如下：画出等高线如下：从两幅图中可以大致看出，当x1=0或1或x2=0这一片区域中，函数值均为0，而点(2，1)处，函数值也为0。下面用不同的初值进行计算首先编制函数的M文件：在命令栏中输入以下命令：输出结果如下：然后改变初值进行计算，输出结果如下：初值最优解函数值迭代次数函数调用次数-1，-1(0.0000，-0.5000)016-2，-2(-0.0003，-1.2880)6.2061e-00722720.1，0.2(0.0000，0.1719)1.1184e-0186243，2(1.8782，1.6809)8.2624e-0119392，0.9(2.0248，0.9066)1.1459e-014530从上表中可以得到以下结论，极小值为0，但是极小点不唯一。当初值选择的不同时，最优解会向x1=0或1或x2=0或者点(2，1)收敛。但是初值不同，收敛的速度不同，初值越接近最优解，收敛越快。下面用不同的搜索方向和步长搜索进行计算在命令栏中输入以下命令(用三种搜索方向(BFGS，DFP和最速下降法)以及两种步长搜索(混合二、三次插值和三次插值)：分析梯度只需将命令中的opt1=optimset(LargeScale,off, MaxFunEvals,1000)改为opt1=optimset(LargeScale,off, MaxFunEvals,1000,GradObj,on)即可。同时函数文件需要作如下改动(此处式子过长不适宜截图)：function 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 nargout1 g(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(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然后在命令栏中输入以下命令：输出结果如下表所示：数值梯度搜索方向步长搜索最优解x最优解y函数值迭代次数调用次数BFGS混合二、三次插值-3.7957e-0071.7187e-0012.9188e-015521DFP-1.9022e-0071.7109e-0017.2776e-016521最速下降1.6709e-0051.6971e-0015.5436e-01218111BFGS三次插值-3.7957e-0071.7187e-0012.9188e-015521DFP-1.9022e-0071.7109e-0017.2776e-016521最速下降1.6709e-0051.6971e-0015.5436e-01218111分析梯度搜索方向步长搜索最优解x最优解y函数值迭代次数调用次数BFGS混合二、三次插值-1.4776e-0051.7567e-0014.5788e-01257DFP-5.7471e-0061.7434e-0016.8437e-01357最速下降1.6238e-0051.7418e-0015.4550e-0122245BFGS三次插值-1.4776e-0051.7567e-0014.5788e-01257DFP-5.7471e-0061.7434e-0016.8437e-01357最速下降1.6238e-0051.7418e-0015.4550e-0122245分析上表，可得到以下结果：1 函数的极小值为0，分别在x1=0或1或x2=0或者点(2，1)处取得。2 当初值不同时，对于得到的最优解、迭代的次数以及函数的调用次数都会有影响。而且当初值不同时，得到的解不同，且得到的都是局部最优解。3 从表中可以看出，使用拟牛顿法的BFGS公式或者DFP公式时，需要的迭代次数相差不多，但是使用最速下降法是需要的迭代系数相对就会较多，函数调用次数也会相应增多。4 分析梯度与数值梯度相比，迭代次数相差不大，相差较多的是函数调用次数。而且在本题中，从数值来看，数值梯度相较于分析梯度要较好一些。下面用自己实现的最速下降法和牛顿法来求解最速下降法在命令栏中输入以下内容：其中前面的p与q是手动求出的梯度，x11=0.5，x21=0.5为初值。输出迭代的结果，x1=0.9614，x2=0.2428，迭代次数为16166次。2 有一组数据ti，yii=1，2，333，其中ti=10(i-1)，yi由下表给出。现要求用这组数据拟合函数fx，t=x1+x2e-x4t+x3e-x5t中的参数x，初值可选为(0.5，1.5，-1，0，0.01，0.02)，用GN和LM两种方法求解。对yi作一扰动，即yi+ei，ei为(-0.05，0.05)内的随机数，观察并分析迭代收敛是否会变慢。iyiiyiiyi10.844120.718230.47820.908130.685240.46730.932140.658250.45740.936150.628260.44850.925160.603270.43860.908170.580280.43170.881180.558290.42480.850190.538300.42090.818200.522310.414100.784210.506320.411110.751220.490330.406初步解决：首先编制函数的M文件：然后在命令栏中输入以下命令：将所有参数输入到MATLAB中。先用LM法计算，输入以下命令：此时，先不使用大规模算法，输出结果如下：然后再在命令栏中输入以下命令：此时使用大规模算法，输出结果如下：从两个输出结果可以看出，当使用了大规模算法之后，迭代次数和函数调用次数要明显少于没有使用大规模算法的时候。但是二者的解一样。再用GN法计算，输入以下命令，同样的，先不使用大规模算法：输出结果如下：然后再在命令栏中输入以下命令：此时使用大规模算法，输出结果如下：从输出结果可以看出，使用了大规模算法之后，迭代次数以及函数调用次数确实减少了很多。但是最终算出的结果是一样的。最终算出的结果如下表所示：x1x2x3x4x50.375411.93585-1.464690.012870.02212迭代次数的比较：不使用大规模算法使用大规模算法LM法358GN法98可以看出，LM法的收敛速度不如GN法快。下面研究当yi发生扰动之后对于迭代收敛速度的影响。(以下全部不使用大规模算法)首先是当ei的所有元素均是相等的时候。ei=-0.4，在命令栏中输入以下内容：LM法：输出结果如下：GN法：输出结果如下：ei=-0.1，在命令栏中输入以下内容：LM法：输出结果如下：GN法：输出结果如下：ei=0.1，在命令栏中输入以下内容：LM法，输出结果如下：GN法，输出结果如下：ei=0.4，在命令栏中输入以下内容：LM法，输出结果如下：GN法，输出结果如下：通过以上结果可以得到下面的表格：结果的比较eix1x2x3x4x50.04LM法0.33541.9358-1.46470.01290.0221GN法0.33541.9358-1.46470.01290.02210.01LM法0.36541.9358-1.46470.01290.0221GN法0.36541.9358-1.46470.01290.02210.01LM法0.38541.9358-1.46470.01290.0221GN法0.38541.9358-1.46470.01290.02210.04LM法0.41541.9358-1.46470.01290.0221GN法0.41541.9358-1.46470.01290.0221迭代速度的比较ei迭代次数0.04LM法35GN法90.01LM法35GN法90.01LM法35GN法90.04LM法35GN法9通过以上的比较可以看出，当yi有一定的小扰动时(ei(-0.05，0.05)，首先发生变化的是x1，从表中可以看出，至始至终在变化的只有x1，而其他的解都没有发生变化。其次，迭代次数与函数调用次数没有发生变化，无论扰动有多大，迭代次数和函数调用次数都还是保持原来的值不变，而且依然是LM法的迭代次数要大于GN法的迭代次数。然后比较ei的元素全部是随机数的情况LM法 在命令栏中输入以下内容：输出结果如下：GN法 在命令栏中输入以下内容：输出结果如下：再换另外一组随机数组进行试验LM法 输出结果如下：GN法 输出结果如下：再换另外一组随机数组进行试验LM法 输出结果如下：GN法 输出结果如下：再换另外一组随机数组进行试验LM法 输出结果如下：GN法 输出结果如下：通过以上结果可以得到下面的表格：结果的比较eix1x2x3x4x51LM法0.30331.0629-0.48490.00780.0287GN法0.30331.0629-0.48500.00780.02872LM法0.34181.1178-0.65460.00980.0342GN法0.34191.1184-0.65530.00980.03413LM法0.2