当前位置:文档之家› 共轭梯度法程序

共轭梯度法程序

一、共轭梯度法共轭梯度法(Conjugate Gradient)是共轭方向法的一种,因为在该方向法中每一个共轭向量都是依靠赖于迭代点处的负梯度而构造出来的,所以称为共轭梯度法。

由于此法最先由Fletcher和Reeves (1964)提出了解非线性最优化问题的,因而又称为FR 共轭梯度法。

由于共轭梯度法不需要矩阵存储,且有较快的收敛速度和二次终止性等优点,现在共轭梯度法已经广泛地应用于实际问题中。

共轭梯度法是一个典型的共轭方向法,它的每一个搜索方向是互相共轭的,而这些搜索方向d仅仅是负梯度方向与上一次迭代的搜索方向的组合,因此,存储量少,计算方便,效果好。

二、共轭梯度法的原理设有目标函数f(X)=1/2X T HX+b T X+c 式1 式中,H作为f(X)的二阶导数矩阵,b为常数矢量,b=[b1,b2,b3,...b n]T 在第k次迭代计算中,从点X(k)出发,沿负梯度方向作一维搜索,得S(K)=-∆f(X(k))式2 X(k+1)=X(k)+ɑ(k)S(k) 式3在式中,ɑ(k)为最优步长。

设与S(k)共轭的下一个方向S(k+1)由点S(k)和点X(k+1)负梯度的线性组合构,即S (k+1)=-∆f (X (k+1))+β(k)S (k) 式4 根据共轭的条件有[S (k)]T ∆2f (X (k))S (k+1)=0 式5 把式2和式4带入式5,得-[∆f(X (k))]T ∆2f (X (k))[-∆f (X (k+1))+β(k)S (k) ]=0 式6 对于式1,则在点X (k)和点X (k+1)的梯度可写为∆f(X (k))=HX (k)+b 式7 ∆f (X (k+1))=HX (k+1)+b 式8 把上面两式相减并将式3代入得ɑ(k)H S (k)=∆f (X (k+1))-∆f(X (k)) 式9 将式4和式9两边分别相乘,并代入式5得-[∆f (X (k+1))+β(k)∆f(X (k))]T [∆f (X (k+1))-∆f(X (k)]=0 式10 将式10展开,并注意到相邻两点梯度间的正交关系,整理后得 β(k )=22||))((||||))1((||k X f k X f ∆+∆ 式11把式11代入式4和式3,得 S (k+1)=-∆f (X (k))+β(k )S (k )X (k+1)=X (k )+ɑ(k )S (k )由上可见,只要利用相邻两点的梯度就可以构造一个共轭方向。

以这种方式产生共轭方向并进行迭代运算的方法,即共轭梯度法。

三、共轭梯度法的迭代方法步骤:1、给定的起始点X (0)和迭代计算精度h ;2、取X (0)的负梯度作为搜索方向S (0)=-∆f (X (0)) 3、沿方向S(k)进行一维搜索,得 X (k+1)=X (k )+ɑ(k )S (k )4、进行收敛检验。

若满足||∆f (X (k+1))||≤h,则令X *=X (k+1),f (X *)=f (X (k+1)) 输出最优解,结束迭代计算;否则,转入(5)5、若k=n ,则令X (0)=X (k+1),转入(2)开始新一轮的迭代;否则,转入(6)。

6、构造新的共轭方向 β=22||))((||||))1((||k X f k X f ∆+∆S (k+1)=-∆f (X (k))+β(k )S (k )令k=k+1,转入(3).四、共轭梯度法的计算框图如下:给定X (0), hS(0)=-∆f (X (0)),k=0minf(X (0))+ɑS (k)) X (k-1)=X (k )+ɑ(k )S (k )||∆f (X(k+1))||≤hk=n ?Β=22||))((||||))1((||k X f k X f ∆+∆S (k+1)=-∆f (X (k))+β(k )S(k )K=k+1X (0)=X (k+1) YNNX *=X (k+1)f (X *)=f (X (k+1))转出Y五、共轭梯度法的C++实现程序清单:#include<stdio.h>#include<math.h>#define N 10#define eps pow(10,-6)double f(double x[],double p[]){double s;s=x[]*x[]+4*p[]*p[]-1;return s;}/*以下是进退法搜索区间源程序*/void sb(double *a,double *b,double x[],double p[]) {double t0,t1,t,h,alpha,f0,f1;int k=0;t0=1; /*初始值*/h=0.01; /*初始步长*/alpha=1;f0=f(x,p,t0);t1=t0+h;f1=f(x,p,t1);while(1){if(f1<f0){h=alpha*h; t=t0; t0=t1; f0=f1; k++;}else{if(k==0){h=-h;t=t1;}else{*a=t<t1?t:t1; *b=t>t1?t:t1; break;}}t1=t0+h;f1=f(x,p,t1);}}double hjfg(double x[],double p[]){double beta,t1,t2,t;double f1,f2;double a=0,b=0;double *c,*d;c=&a,d=&b;sb(c,d,x,p);/*调用进退法搜索区间*/printf("\nx1=%lf,x2=%lf,p1=%lf,p2=%lf",x[0],x[1],p[0],p[1]); printf("\n[a,b]=[%lf,%lf]",a,b);beta=(sqrt(5)-1.0)/2;t2=a+beta*(b-a); f2=f(x,p,t2);t1=a+b-t2; f1=f(x,p,t1);while(1){if(fabs(t1-t2)<eps)break;else{if(f1<f2){t=(t1+t2)/2;b=t2; t2=t1;f2=f1; t1=a+b-t2; f1=f(x,p,t1);}else{a=t1; t1=t2;f1=f2;t2=a+beta*(b-a);f2=f(x,p,t2);}}}t=(t1+t2)/2;return t;}/*以下是共轭梯度法程序源代码*/ void gtd(){double x[N],g[N],p[N],t=0,f0,mod1=0,mod2=0,nanda=0;int i,k,n;printf("请输入函数的初值x1,x2=");scanf("%d",&n);printf("\n请输入步长值:\n");for(i=0;i<n;i++)scanf("%lf",&x[i]);f0=f(x,g,t);g[0]=2*x[0]; g[1]=50*x[1];mod1=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的长度*/if(mod1>eps){p[0]=-g[0]; p[1]=-g[1]; k=0;while(1){t=hjfg(x,p);printf("\np1=%lf,p2=%lf,t=%lf",p[0],p[1],t);x[0]=x[0]+t*p[0]; x[1]=x[1]+t*p[1];g[0]=2*x[0]; g[1]=50*x[1];/*printf("\nx1=%lf,x2=%lf,g1=%lf,g2=%lf",x[0],x[1] ,g[0],g[1]);*/mod2=sqrt(pow(g[0],2)+pow(g[1],2)); /*求梯度的长度*/if(mod2<=eps) break;else{if(k+1==n){g[0]=2*x[0]; g[1]=50*x[1];p[0]=-g[0]; p[1]=-g[1]; k=0;}else{nanda=pow(mod2,2)/pow(mod1,2);printf("\nnanda=%lf,mod=%lf",nanda,mod2);p[0]=-g[0]+nanda*p[0];p[1]=-g[1]+nanda*p[1];mod1=mod2;k++;}}printf("\n--------------------------");}}printf("\n最优解为x1=%lf,x2=%lf",x[0],x[1]);printf("\n最终的函数值为%lf",f(x,g,t));}main(){gtd();}五、共轭梯度法的函数:用共轭梯度法求f(X)=x21+4x22-1的极小值,设X(0)=[1,1]T,h=0.01.六、函数运行:运行结果如下:请输入初始值:1 1请输入步长值:0.01最优解为x1=-0.000000,x2=-0.000000最终的函数值为-1Press any key to continue七、小结经过此次优化设计学习及作业完成,我基本掌握了优化设计课程的知识内容及任务,特别是在老师安排下动手完成了共轭梯度法的编程函数实现,对 C++语言也更熟悉,对以后的进一步的学习有很大的帮助,尽管在编程时出现诸多错误,甚至无法实现函数功能,但在老师不断地帮助下及时改进,终于实现函数,第一次将C++用于解决实际问题,收获很大。

相关主题