运筹学实验共轭梯度法.doc

共轭梯度法一、实验目的(1).熟悉使用共轭梯度法求解无约束非线性规划问题的原理；(2).在掌握原理的基础上熟练运用此方法解决问题；(3).学会利用计算机语言编写程序来辅助解决数学问题；(4).解决问题的同时分析问题，力求达到理论与实践的相统一；(5).编写规范的实验报告.二、问题描述minfx=x12+x22-x1x2+2x1-4x2自选初始点开始迭代三、算法介绍：共轭梯度法为求解线性方程组而提出。后来，人们把这种方法用于求解无约束最优化问题，使之成为一种重要的最优化方法。共轭梯度法的基本思想是把共轭性与最速下降方法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜索，求出目标函数的极小点。根据共轭方向的基本性质，这种方法具有二次终止性。在各种优化算法中，共轭梯度法是非常重要的一种。其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。共轭方向无约束最优化方法的核心问题是选择搜索方向.在本次实验中,我们运用基于共轭方向的一种算法共轭梯度法：给定初始点(0,0)，k=1，最大迭代次数n确定搜索方向进退法确定搜索区间分割法确定最优步长四、程序%,function m,k,d,a,X,g1,fv = GETD( G,b,c,X,e,method)if nargin=e k(i-1)=(m/m1)2; d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); %a1(i)=-(X(:,i)*G*d(:,i)+b*d(:,i)/(d(:,i)*G*d(:,i);a(i)=g1(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=g1 subs(subs(g,x1,X(1,i+1),x2,X(2,i+1); m1=m; m=norm(g1(:,i+1); i=i+1; end case PRP while m=e k(i-1)=g1(:,i)*(g1(:,i)-g1(:,i-1)/(norm(g1(:,i-1)2; d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=g1 subs(subs(g,x1,X(1,i+1),x2,X(2,i+1); m=norm(g1(:,i+1); i=i+1; end case HS while m=e k(i-1)=g1(:,i)*(g1(:,i)-g1(:,i-1)/(d(:,i-1)*(g1(:,i)-g1(:,i-1); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=g1 subs(subs(g,x1,X(1,i+1),x2,X(2,i+1); m=norm(g1(:,i+1); i=i+1; end case DY while m=e k(i-1)=g1(:,i)*g1(:,i)/(d(:,i-1)*(g1(:,i)-g1(:,i-1); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=g1 subs(subs(g,x1,X(1,i+1),x2,X(2,i+1); m=norm(g1(:,i+1); i=i+1; end case LS while m=e k(i-1)=g1(:,i)*(g1(:,i)-g1(:,i-1)/(d(:,i-1)*(-g1(:,i-1); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); %a(i)=-(X(:,i)*G*d(:,i)+b*d(:,i)/(d(:,i)*G*d(:,i); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=g1 subs(subs(g,x1,X(1,i+1),x2,X(2,i+1); m=norm(g1(:,i+1); i=i+1; end case CD while m=e k(i-1)=g1(:,i)*g1(:,i)/(d(:,i-1)*(-g1(:,i-1); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=g1 subs(subs(g,x1,X(1,i+1),x2,X(2,i+1); m=norm(g1(:,i+1); i=i+1; end case WYL while m=e k(i-1)=g1(:,i)*(g1(:,i)-(m/m1)*g1(:,i-1)/(m12); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)*g1(:,i)/(d(:,i)*G*d(:,i); %a(i)=-(X(:,i)*G*d(:,i)+b*d(:,i)/(d(:,i)*G*d(:,i); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=g1 subs(subs(g,x1,X(1,i+1),x2,X(2,i+1); m1=m; m=norm(g1(:,i+1); i=i+1; endendfv=subs(subs(f,x1,X(1,i),x2,X(2,i); endl1=X(1,i);l2=X(2,i);w1=X(1,1);w2=X(2,1);v1=min(l1,w1)-abs(l1-w1)/10:abs(l1-w1)/10:max(l1,w1)+abs(l1-w1)/10;v2=min(l2,w2)-abs(l2-w2)/10:abs(l2-w2)/10:max(l2,w2)+abs(l2-w2)/10;x,y=meshgrid(v1,v2);s=size(x);z=zeros(size(x);for i=1:s(1) for j=1:s(2) z(i,j)=1/2*x(i,j),y(i,j)*G*x(i,j);y(i,j)+b*x(i,j);y(i,j)+c; endendpx,py = gradient(z,.2,.2);contour(v1,v2,z), hold on, quiver(v1,v2,px,py)C,h = contour(x,y,z);set(h,ShowText,on,TextStep,get(h,LevelStep)*2)x1=X(1,:);y1=X(2,:);plot(x1,y1,r*:); 五、计算结果如图所示，输入 G=2,-1;-1,2; b=2;-4; c=0;X=1;1;e=1e-3;method=FR;表示正定矩阵为G，b为一次元系数矩阵，c为常数0，X是初始迭代点，e为精度，method为共轭梯度的方法。得到的结果如图为：迭代一次就出来结果了0;2,同时此时函数的最小值为-4同时得到共轭梯度收敛时候的图：六、结果分析由FR法得到的最优解为（0，2），迭代的次数为1，由于精度的降低（= ），实际证明对结果也有一定的影响。在本次实验中取得的最优值为-4，但可能因为初值点取得比较特殊，如果换成1；0.5时：结果虽然无限接近0;2，迭代的次数也变为2，所以由此可以得出一个好的初值点能加快现实中的求解，不过毕竟这是电脑的工作，所以一笑而过。七、心得体会本题不是很难，根据算法过程编写程序需要的是耐心和仔细，本问题要考虑的方面比较多，