最优化数值实验报告一

最优化数值实验报告一 实验目的: 1，掌握对无约束最优化问题的各种解法，主要有最速下降法，牛顿法，共轭梯度法，变尺度法(DFP和BFGS法)，非线性最小二乘法。 2，用MATLAB变成求解问题，并对各种方法的计算过程和结果进行分析。 实验原理:(略) 实验内容: (结合了《最优化练习题017》中的第2题进行分析。) 求解下面的方程组 ,xxx3,cos(),1/2,0,123,2fxxxx(),,81(，0.1)，sin，1.06,0 ,123 ,1,xx12ex,，20，(10,3),03,3, 问题分析: 对于上述问题，可以用两大类的方法求解。其一是将上述方程组转化为无约束的最小二乘问题，并非别用最速下降法，牛顿法，共轭梯度法，变尺度法，非线性最小二乘法求解。其二是将上述方程组转化为有约束的问题，并用惩罚函数法和广义乘子法进行求解。下面我们将上述问题转化为无约束的问题进行求解。 将上述问题转化为最小二乘问题即是对上述问题进行转化，而变成平方和的形式，我们 minf(x)f(x)可以将原问题转为求的无约束优化问题，其中如下， 1,xx222212f(x),(3x,cos(xx),1/2)，(x,81(x，0.1)，sinx，1.06)，(e，20x，(10,,3))12312333 下面我们借助MATLAB用不同的优化方法进行分析求解。 另外，在下面的的搜索算法中，当需要做一维搜索时，都是用进退法和黄金分割法来实现的，其中进退法确定一个单峰区间，而黄金分割法则进行精确的一维搜索。(分别对应程序中的jintuifa.m和gold.m,具体程序见附录)。值得指出的是，为了保证能够使黄金分割法正常工作，进退法中步长选取显得很是重要，因为黄金分割法主要是针对单峰区间来进行搜索的。步长选择的太大，难以确定单峰区间;步长选择的太小，计算量就比较大，迭代的较慢。下面实验中选择的步长为0.0125。 1,最速下降法 最速下降法，通俗点来说，就是选取一个目标函数值下降最快的方法，以利于尽快地达到极小点。最速下降法的关键就是最速下降方向的选取，实际上我们知道负梯度方向为最速下降方向。然后我们进行一维搜索，直到满足精度要求，则停止计算。 d,,实验中，我们一律取初始点为(0，0，0)，当搜索方向d满足时，停止搜索， ,6,其中= 。 1*10 程序: syms x1 x2 x3 r; f=(3*x1-cos(x2*x3)-1/2)^2+(x1^2-81*(x2+0.1)+sin(x3)+1.06)^2+(exp(-x1*x2)+20 *x3+1/3*(10*3.14159-3))^2; df=(jacobian(f,[x1 x2 x3])).';%对f求导 x0=[0 0 0]';error=1e-6; %取初始点为(0，0，0);允许误差为1e-6 g=subs(df,[x1 x2 x3],x0.'); g=double(g);k=0; while(norm(g)>error) d=-g; %确定最优搜索方向 z=subs(f,[x1 x2 x3],(x0+r*d).'); result=jintuifa(z,r); %进退法确定搜索区间 result2=gold(z,r,result); %黄金分割法确定最优步长 step=result2(1); x0=x0+step*d; g=subs(df,{x1 x2 x3},{x0(1,1),x0(2,1),x0(3,1)}); g=double(g); k=k+1; end; k x0 min=subs(f,[x1 x2 x3],x0) 结果: 由最速下降法求得的最优解为(0.4996，-0.0900，-0.5259)，并且迭代次数为2129，最优值2.4829e-014，显然迭代的次数太多，花费的时间太多，所以效率并不高。 究其原因，应该是锯齿现象的影响。也就是说，虽然从局部看，最速下降算法是函数值下降最快的方向，但从全局看，由于锯齿现象，使得向极点移动的过程中要总弯路，因此收敛速度变慢。 2，牛顿法 牛顿法是通过x(k)和在这一点的目标函数的梯度和HESSE矩阵的逆来得到后续点x(k+1),在适当的条件下，这个序列x(k)将收敛。 kk，1,6,x,x,, 同样，在实验中我们取初始点(0，0，0)，当，停止寻优，其中= 。 1*10 程序: syms x1 x2 x3; f=(3*x1-cos(x2*x3)-1/2)^2+(x1^2-81*(x2+0.1)+sin(x3)+1.06)^2+(exp(-x1*x2)+20* x3+1/3*(10*3.14159-3))^2; x=[x1 x2 x3]; x0=[0 0 0]';error=1e-6;t=1;m=0; while (norm(t)>error) a=subs(diff(f,x1),x,x0);b=subs(diff(f,x2),x,x0);c=subs(diff(f,x3),x,x0); df=[a;b;c]; d=subs(diff(f,x1,2),x,x0);h=subs(diff(f,x2,2),x,x0); l=subs(diff(f,x3,2),x,x0);e=subs((diff(diff(f,x1),x2)),x,x0); f1=subs((diff(diff(f,x1),x3)),x,x0);g=subs((diff(diff(f,x2),x1)),x,x0); i=subs(diff(diff(f,x2),x3),x,x0);j=subs(diff(diff(f,x3),x1),x,x0); k=subs(diff(diff(f,x3),x2),x,x0); df2=[d e f1;g h i;j k l];%HESSE矩阵 df=double(df);df2=double(df2); if(det(df2)==0) break; else xx=x0-inv(df2)*df; end t=xx-x0;x0=xx;m=m+1; end m x0 min=subs(f,[x1 x2 x3],x0) 结果: 由牛顿法求得的最优解为(0.4996，-0.0900，-0.5259)，并且迭代次数为6，显然牛顿法的收敛速度很快，最优值为8.0119e-031。但是牛顿法由于要求HESSE矩阵的逆，而计算中并不能保证HESSE矩阵总是可逆的，所以，当初始点选取一定值时，就会出现问题，比如在本次实验中取初始点为(100，100，100)，就出现了问题，这也是牛顿法最主要的缺陷。同时它也不能保证每次的搜索方向为下降方向。 3，共轭梯度法 无约束问题的核心就是搜索方向的选取，共轭梯度法的基本思想就是把共轭性和最速下降法想结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向就行搜索，求出目标函数最小点。实验中采用的是FR法来构造共轭方向。 k，1,f(x),, 在实验中取初始点(0，0，0)，当时，停止寻优，由于计算量偏大，这 ,6,次实验中选取的= 。 1*10 程序: syms x1 x2 x3 r; f=(3*x1-cos(x2*x3)-1/2)^2+(x1^2-81*(x2+0.1)+sin(x3)+1.06)^2+(exp(-x1*x2)+20* x3+1/3*(10*3.14159-3))^2; x=[x1,x2,x3];df=jacobian(f,x);df=df.'; error=0.000001;x0=[0,0,0]';g1=subs(df,x,x0);k=0; while(norm(g1)>error) if k==0 d=-g1; else bta=g1'*g1/(g0'*g0); d=-g1+bta*d0; end y=subs(f,x,x0+r*d); result=jintuifa(y,r); result2=gold(y,r,result); step=result2; x0=x0+step*d; g0=g1;g1=subs(df,x,x0); d0=d;k=k+1; end; k x0 min=subs(f,[x1 x2 x3],x0) 结果: 由FR法得到的最优解为(0.4996, -0.0900, -0.5259)，迭代的次数为32，由于精度的 ,6,降低(= )，实际 证明 住所证明下载场所使用证明下载诊断证明下载住所证明下载爱问住所证明下载爱问 对结果也有一定的影响。在本次实验中取得的最优值为1*10 1.1496e-018，而前面牛顿法迭代6次算出的最优值为8.0119e-031，效果要好了很多。所以说，收敛速度相对比较慢，应该是它的一个缺点。 4，DFP法 牛顿法的收敛的速度很快，前面我们已经发现了，但是这都要保证HESSE可逆才可以，所以就有拟牛顿法，它不用求逆，而用其它的矩阵来近似。拟牛顿法是一种比较有效的算法。DFP也提出了一种这样的近似矩阵。具体的构造方法程序中可以看到。 k，1,6,,f(x),, 同样，在实验中我们取初始点(0，0，0)，当，停止寻优，其中= 。 1*10 程序: syms x1 x2 x3 r; f=(3*x1-cos(x2*x3)-1/2)^2+(x1^2-81*(x2+0.1)+sin(x3)+1.06)^2+(exp(-x1*x2)+20* x3+1/3*(10*3.14159-3))^2; x=[x1,x2,x3];df=jacobian(f,x);df=df.'; error=1e-6;x0=[0,0,0]';g1=subs(df,x,x0);k=0;H0=[1,0 0;0,1 0;0 0 1]; while(norm(g1)>error) if k==0 d=-H0*g1; else H1=H0+pk*pk'/(qk'*pk)-H0*qk*qk'*H0/(qk'*H0*qk); d=-H1*g1;H0=H1; end z=subs(f,x,x0+r*d); result=jintuifa(z,r); result2=gold(z,r,result); step=result2; x0=x0+step*d; g0=g1;g1=subs(df,x,x0); qk=g1-g0; pk=step*d; k=k+1; end; k x0 min=subs(f,[x1 x2 x3],x0) 结果: 由DFP法得到的最优解为(0.4996，-0.0900，-0.5259)，迭代的次数为7，速度还是非常快的。取得的最优值为4.2298e-019。只要精度取得足够高，DFP还是很理想的。但是，计算机计算过程中的舍入误差的影响也是不可以忽视的，因为这会导致校正矩阵计算过程中出现除数为0的情况，实际上在实验的过程中也碰到了这样的问题。 5，BFGS法 所谓BFGS法，也是变尺度法的一种，对DFP算法做了一点修正。 程序: syms x1 x2 x3 r; f=(3*x1-cos(x2*x3)-1/2)^2+(x1^2-81*(x2+0.1)+sin(x3)+1.06)^2+(exp(-x1*x2)+20* x3+1/3*(10*3.14159-3))^2; x=[x1,x2,x3];df=jacobian(f,x);df=df.'; error=1e-6;x0=[0,0,0]';g1=subs(df,x,x0);k=0;H0=[1,0 0;0,1 0;0 0 1]; while(norm(g1)>error) if k==0 d=-H0*g1; else H1=H0+(1+qk'*H0*qk/(pk'*qk))*(pk*pk')/(pk'*qk)-(pk*qk'*H0+H0*qk* pk')/(pk'*qk); d=-H1*g1;H0=H1; end z=subs(f,x,x0+r*d); result=jintuifa(z,r); result2=gold(z,r,result); step=result2; x0=x0+step*d; g0=g1;g1=subs(df,x,x0); qk=g1-g0; pk=step*d; k=k+1; end; k x0 min=subs(f,[x1 x2 x3],x0) 结果: 由BFGS法得到的最优解为(0.4996，-0.0900，-0.5259)，迭代的次数为7，速度还是非常快的。取得的最优值为9.4951e-020。从结果来看，它比DFP法略微要好。 6，非线性最小最小二乘法 程序: syms x1 x2 x3 r; f=(3*x1-cos(x2*x3)-1/2)^2+(x1^2-81*(x2+0.1)+sin(x3)+1.06)^2+(exp(-x1*x2)+20 *x3+1/3*(10*3.14159-3))^2; f1=3*x1-cos(x2*x3)-1/2;f2=x1^2-81*(x2+0.1)+sin(x3)+1.06;f3=exp(-x1*x2)+20*x 3+1/3*(10*3.14159-3); y=[f1;f2;f3];x=[x1,x2,x3]; df=jacobian([f1;f2;f3],x); df=df.'; error=1e-6;x0=[0,0,0]';Ak=subs(df,x,x0);fk=subs(y,x,x0);k=0;t=[1 1 1]';j=0; while(norm(t)>error) Bk=Ak'*Ak; if(det(Bk)==0) break; else d=-1*inv(Bk)*Ak'*fk; end z=subs(f,x,x0+r*d); result=jintuifa(z,r); result2=gold(z,r,result); step=result2; x0=x0+step*d; t=step*d; k=k+1; end; k x0 min=subs(f,[x1 x2 x3],x0) 结果: 由非线性最小二方法得到的最优解为(0.5098，-0.0886，-0.5294)，迭代的次数为2， TAA速度是最快的。取得的最优值为0.0174。但是，它也有着局限性。并不能保证()总 是可逆的。 实验结论: 上面我们已经对几种不同的优化方法进行了比较分析。在每一部分已经对各种方法的优 缺点进行了分析。在具体的实际应用中，要根据不同的情况选择合适的方法。 附录 进退法确定一维搜索区间 function result=jintuifa(y,r) t0=0; step=0.0125; t1=t0+step; ft0=subs(y,{r},{t0});ft1=subs(y,{r},{t1}); if(ft1<=ft0) step=2*step;t2=t1+step;ft2=subs(y,{r},{t2}); while(ft1>ft2) t1=t2;step=2*step;t2=t1+step;ft1=subs(y,{r},{t1});ft2=subs(y,{r},{t2}); end else step=step/2;t2=t1;t1=t2-step;ft1=subs(y,{r},{t1}); while(ft1>ft0) step=step/2;t2=t1;t1=t2-step;ft1=subs(y,{r},{t1}); end end result=[t2]; 黄金分割法进行一维搜索 function result=gold(y,r,m) a=0;b=m;e=1e-5;a1=a+0.382*(b-a);f1=subs(y,{r},{a1});a2=a+0.618*(b-a);f2=subs(y,{r},{a2}); while abs(b-a)>=e if f1