Conjugate-gradient-methodPBFGS.rar_BFGS算法C语言_梯度优化

/* * BFGS.c * * Created on: 2009-6-4 * Author: Administrator */ //×××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××× //××××××××××××××××××××××××××××BFGS算法C程序××××××××××××××××××××××××××××××× //×××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××× //×××××××××××本程序适用于2设计变量的函数优化问题，对于不同的设计变量××××××××× //***********个数可以改变维数,对于DFP算法只需修改校正矩阵即可。对于×××××××××× //×××××××××××BFGS算法的C++程序与之类似。 ×××××××× //*××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××× //×××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××× //************************************************************************ #include<stdio.h> #include<stdlib.h> #include<math.h> #include<conio.h> #define tt 0.01 //-------------------一维搜索初始步长--------------------------- #define ff 1.0e-6 //-------------------差分法求梯度时的步长----------------------- #define ac 1.0e-6 //-------------------终止迭代梯度模收敛精度--------------------- #define ad 1.0e-6 //-------------------一维搜索收敛精度--------------------------- #define ax 1.0e-6 //-------------------终止迭代点距收敛精度----------------------- #define n 2 //-------------------设计变量的维数----------------------------- double ia; /*FILE *fp;*/ //=================================输入目标函数表达式======================================== double fny(double *x) { double x1=x[0],x2=x[1]; double f; f=x1*x1+2*x2*x2-4*x1-2*x1*x2; return f; } //========================================================================================== //==================================迭代更新最优点=========================================== double * iterate(double *x,double a,double *s) { double *x1; int i; x1=(double *)malloc(n*sizeof(double)); for(i=0;i<n;i++) x1[i]=x[i]+a*s[i]; return x1; } //========================================================================================== //===================================计算更新后节点函数值==================================== double func(double *x,double a,double *s) { double *x1; double f; x1=iterate(x,a,s); f=fny(x1); return f; } //========================================================================================== //===================================确定初始单谷区间的退进法================================ void finding(double a[3],double f[3],double *xk,double *s) { double t=tt; int i; double a1,f1; a[0]=0;f[0]=func(xk,a[0],s); for(i=0;;i++) //循环次数i，受break执行的影响 { a[1]=a[0]+t; f[1]=func(xk,a[1],s); if(f[1]<f[0]) break; if(fabs(f[1]-f[0])>=ad) { t=-t; a[0]=a[1];f[0]=f[1]; } else { if(ia==1) return; //break t=t/2;ia=1; } } for(i=0;;i++) { a[2]=a[1]+t; f[2]=func(xk,a[2],s); if(f[2]>f[1]) break; t=2*t; a[0]=a[1];f[0]=f[1]; a[1]=a[2];f[1]=f[2]; } if(a[0]>a[2]) { a1=a[0]; f1=f[0]; a[0]=a[2]; f[0]=f[2]; a[2]=a1; f[2]=f1; } return; } //====================================================================================== //=======================================一维搜索======================================== double lagrange(double *xk,double *ft,double *s) { int i; double a[3],f[3]; double b,c,d,aa; finding(a,f,xk,s); for(i=0;;i++) { if(ia==1) { aa=a[1]; *ft=f[1]; break; } d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2]); if(fabs(d)==0) break; c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d; if(fabs(c)==0) break; b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1]); aa=-b/(2*c); *ft=func(xk,aa,s); if(fabs(aa-a[1])<=ad) {if(*ft>f[1]) aa=a[1];break;} if(aa>a[1]) { if(*ft>f[1]) {a[2]=aa;f[2]=*ft;} else if(*ft<f[1]) {a[0]=a[1];a[1]=aa;f[0]=f[1];f[1]=*ft;} else if(*ft==f[1]) { a[2]=aa;a[0]=a[1]; f[2]=*ft;f[0]=f[1]; a[1]=(a[0]+a[2])/2; f[1]=func(xk,a[1],s); } } else { if(*ft>f[1]) {a[0]=aa;f[0]=*ft;} else if(*ft<f[1]) {a[2]=a[1];a[1]=aa;f[2]=f[1];f[1]=*ft;} else if(*ft==f[1]) {a[0]=aa;a[2]=a[1]; f[0]=*ft;f[2]=f[1]; a[1]=(a[0]+a[2])/2; f[1]=func(xk,a[1],s); } } } if(*ft>f[1]) {*ft=f[1];aa=a[1];} return aa; } //======================================================================================= //====================================计算向量梯度======================================== double *gradient(double *xk) { double *g,f1,f2,q; int i; g=(double*)malloc(n*sizeof(double)); f1=fny(xk); for(i=0;i<n;i++) {q=ff; xk[i]=xk[i]+q; f2=fny(xk); g[i]=(f2-f1)/q; xk[i]=xk[i]-q; } return g; } //======================================================================================== //=======================================BFGS算法主程序======================================================= double * bfgs(double *xk) { double u[n],v[n],h[n][n],dx[n],dg[n],s[n]; double aa,ib; double *ft,*xk1,*g1,*g2,*xx,*x0=xk; double fi,fxk1; int i,j,k,step=1; ft=(double *)malloc(sizeof(double)); xk1=(double *)malloc(n*sizeof(double)); //=====================================数据初始化============================================== for(i=0;i<n;i++) //--------------尺度矩阵h[n][n]、搜索方向s[n]初始化-------- { s[i]=0; for(j=0;j<n;j++) { h[i][j]=0; if(j==i) h[i][j]=1; } } //--------------------------------------------------------- g1=gradient(xk); //----------------向量梯度初始化g1[n]---------------------- fi=fny(xk); x0=xk; printf("%d %f %f %f \n",step,x0[0],x0[1],fi); //============================================================================================ //=========================================迭代求解最优点====================================== for(k=0;k<n;k++) { ib=0; if(ia==1) { xx=xk; break; } ib=0; for(i=0;i<n;i++) s[i]=0; //----------------求搜索方向s[n]--------------------- for(i=0;i<n;i++) for(j=0;j<n;j++) s[i]+= -h[i][j]*g1[j]; //--------------------------------------------------- aa=lagrange(xk,ft,s); //----------------计算步长a-------------------------- xk1=iterate(xk,aa,s); //----------------迭代计算X[n]----------------------- g2=gradient(xk1); //----------------计算向量梯度g[n]------------------- fxk1=fny(xk1); step+=1; printf("%d %f %f %f \n",step,xk1[0],xk1[1],fxk1); if(sqrt((g2[1]-g1[1])*(g2[1]-g1[1])+(g2[0]-g1[0])*(g2[0]-g1[0]))>=ac||sqrt((xk1[1]-xk[1])*(xk1[1]-xk[1])+(xk1[0]-xk[0])*(xk1[0]-xk[0]))>=ax) //---------------------判断迭代是否结束---------------------- {ib=ib+1;} if(ib==0) { xx=xk1; break; } //---------------------------------------- fi=*ft; if(k==0)//------------------------------求解校正矩阵------------------------------- { int j; double a1=0,a2=0; for(i=0;i<n;i++) { dg[i]=g2[i]-g1[i]; dx[i]=xk1[i]-xk[i]; } for(i=0;i<n;i++) { int j; u[i]=0;v[i]=0; for(j=0;j<n;j++) { u[i]=u[i]+dg[j]*h[j][i]; v[i]=v[i]+dg[j]*h[i][j]; } } for(j=