dfp_bfgs.cpp

无约束优化DFP&BFGS方法的C++实现

#include <iostream.h>
#include <math.h>
/*#define	kethe 1e-5*/
#define	n 2

//=========分配双精度二维动态数组==========================================
void Make2DArray_double(double **&P, int ne, int nL)
{
    P=new double *[ne];
	for(int i=0;i<ne;i++)P[i]=new double[nL];
}
//=========分配双精度二维动态数组==========================================
//=========用于一维搜索的原函数============================================
double f(double *x0,double *S,double alphai)
{//POWELL法的原函数
	double x[n],f1;
	for(int i=0;i<n;i++)x[i]=x0[i]+alphai*S[i];
//	f1=pow(x[0],2)+2*pow(x[1],2)-4*x[0]-2*x[0]*x[1];
//	f1=10*pow(x[0]+x[1]-5,2)+pow(x[0]-x[1],2);
//	f1=pow(x[0]-5,2)+pow(x[1]+10,2)+pow(x[2]+3,2)+x[0]*x[2]+7;
	f1=2*pow(x[0],2)+3*pow(x[1],2)-8*x[0]+10;
//	f1=4*pow(x[0]-5,2)+pow(x[1]-6,2);
	return(f1);
}
//=========原函数==========================================================
//=========用于一维搜索的原函数============================================
double f2(double *x)
{//原函数
	double f1;
//	f1=pow(x[0],2)+2*pow(x[1],2)-4*x[0]-2*x[0]*x[1];
//	f1=10*pow(x[0]+x[1]-5,2)+pow(x[0]-x[1],2);
//	f1=pow(x[0]-5,2)+pow(x[1]+10,2)+pow(x[2]+3,2)+x[0]*x[2]+7;
	f1=2*pow(x[0],2)+3*pow(x[1],2)-8*x[0]+10;
//	f1=4*pow(x[0]-5,2)+pow(x[1]-6,2);
	return(f1);
}
//=========原函数==========================================================
//=========计算梯度========================================================
void calculate_g(double *x, double *gk)
{
//	gk[0]=8*(x[0]-5);
//	gk[1]=2*(x[1]-6);
	gk[0]=4*x[0]-8;
	gk[1]=6*x[1];
}
//=========计算梯度========================================================
//=========计算范数========================================================
double fanshu(double *gk)
{
	int i;
	double fanshu=0;
	for (i=0;i<n;i++)fanshu=fanshu+pow(gk[i],2);
	fanshu=sqrt(fanshu);
	return fanshu;
}
//=========计算范数========================================================
//=========初始化==========================================================
void Initialize(double *Xk, double *kethe, double *alpha0, double *h,double *L)
{
	int i;
	for (i=0;i<n;i++){
		Xk[i]=i+8;
	}
	*kethe=1e-5;
	*alpha0=1.5;
	*h=0.1;
	*L=1e-3;
}
//=========初始化==========================================================
//=========初始化Ak========================================================
void Initialize_Ak(double **Ak)
{
	int i,j;
	for (i=0;i<n;i++) {
		for (j=0;j<n;j++){
			if(i==j)Ak[i][j]=1;
			else Ak[i][j]=0;	
		}
	}
}
//=========初始化Ak========================================================
//=========排列大小========================================================
void arrange(double *x1,double *x3)
{
	double temp;
	if(*x1>*x3)
	{
		temp=*x1;
		*x1=*x3;
		*x3=temp;
	}
}
//=========排列大小========================================================
//=========进退法==========================================================
void jintui(double x0, double h, double *x1,double *x2,double *x3,double *x,double *S)
{
	double x4;
	double f1,f2,f4;
	int k=0;
	*x1=0;*x2=0;*x3=0;x4=0;
	f1=0;f2=0;f4=0;
	*x1=x0;
	f1=f(x,S,*x1);
	for(k=1; ;k++)
	{
		x4=*x1+h;
		f4=f(x,S,x4);
		if(f(x,S,x4)<f(x,S,*x1))
		{//前进
			*x2=*x1;
			*x1=x4;
			h=2*h;
			x4=*x1+h;
			if((f(x,S,*x2)-f(x,S,*x1))*(f(x,S,*x1)-f(x,S,x4))<=0)
			{//进退终止条件
				*x3=*x2;
				*x2=*x1;
				*x1=x4;
				arrange(x1,x3);
				break;
			}
			continue;
		}
		else
		{//后退
			if(k==1)
			{
				h=-h;
				*x2=x4;
				f2=f(x,S,x4);
				continue;
			}
			else
			{				
				if((f(x,S,*x2)-f(x,S,*x1))*(f(x,S,*x1)-f(x,S,x4))<=0)
				{//进退终止条件
					*x3=*x2;
					*x2=*x1;
					*x1=x4;
					arrange(x1,x3);
					break;
				}
			}
		}	
	}
}
//========进退法===========================================================
//========黄金分割法=======================================================
double f0618f(double a1, double b1, double L,double *x,double *S)
{
	double a2,b2,lamda1,lamda2,mu1,mu2,h;
	int k=1;
	h=(pow(5,0.5)-1)/2;
	lamda1=a1+(1-h)*(b1-a1);	//初始化试探点
	mu1=a1+h*(b1-a1);			//初始化试探点	
	for(k=0;b1-a1>=L ;k++)
	{
		if(f(x,S,lamda1)>f(x,S,mu1))
		{
			a2=lamda1;
			b2=b1;
			lamda2=mu1;
			mu2=a2+h*(b2-a2);
		}
		else
		{
			a2=a1;
			b2=mu1;
			mu2=lamda1;
			lamda2=a2+(1-h)*(b2-a2);
		}
		a1=a2;
		b1=b2;
		lamda1=lamda2;
		mu1=mu2;
	}
	return((a1+b1)/2);
}
//========黄金分割法=======================================================
void main()
{
	int i,j,k,t,w=0;
	double *Xk=new double [n];		//用于存储变量计算点
	double *Xk1=new double [n];		//用于存储修正后的变量计算点
	double *gk=new double [n];		//用于存储梯度
	double *gk1=new double [n];		//用于存储下一个迭代点的梯度
	double *Sk=new double [n];		//用于存储搜索方向
	double *sigma=new double [n];	//用于存储sigma
	double *yk=new double [n];		//用于存储yk
	
//	double *alpha_k=new double [n];	//用于存储最佳步长
	double **Ak, **Ak1;				//用于存储构造矩阵
	double **Mk, **Nk;
	double **temp7, **temp8;
	Make2DArray_double(Ak,n,n);
	Make2DArray_double(Mk,n,n);
	Make2DArray_double(Nk,n,n);
	Make2DArray_double(temp7,n,n);
	Make2DArray_double(temp8,n,n);
	double kethe,alpha0,h,L;		//用于存储精度和一维搜索的初始点以及步长
	double x1,x2,x3,temp1;
	double *temp4=new double [n];
	double *temp5=new double [n];
	double *temp6=new double [n];
	double alpha_k;					//用于存储最佳步长

	Initialize(Xk, &kethe, &alpha0,&h,&L);				//初始化
	calculate_g(Xk, gk);								//计算梯度
	for(w=0; ;w++)
	{
		Initialize_Ak(Ak);								//初始化Ak
		for(k=0;k<n;k++) 
		{
			cout<<"w="<<w<<"		k="<<k<<endl;
			cout<<"搜索方向: ";
			for (i=0;i<n;i++) {							//计算搜索方向
				Sk[i]=0;//初始化Sk
				for (j=0;j<n;j++) Sk[i]=Sk[i]+Ak[i][j]*gk[j];
				Sk[i]=(-1)*Sk[i];
				cout<<"Sk[i]="<<Sk[i]<<"	";
			}cout<<endl;
			
			jintui(alpha0,h,&x1,&x2,&x3,Xk,Sk);			//一维搜索求各个方向上的最佳步长所在的单峰区间
			alpha_k=f0618f(x1,x3,L,Xk,Sk);				//一维搜索求各个方向上的最佳步长
			cout<<"最佳步长alpha_k为："<<alpha_k<<endl;
			for(i=0;i<n;i++)Xk1[i]=Xk[i]+alpha_k*Sk[i];	//计算下一个迭代点



			cout<<"迭代点: ";
			for (i=0;i<n;i++) cout<<"Xk1[i]="<<Xk1[i]<<"		";	
			cout<<endl;

			calculate_g(Xk1, gk1);						//计算下一个迭代点的梯度


			if (fanshu(gk1)<kethe)						//如果满足精度条件，则结束迭代
			{
				cout<<endl<<"fanshu(gk1)="<<fanshu(gk1)<<"	w="<<w<<endl;
				w=-1;
				break;				
			}
			else
			{											//否则用DFP公式计算构造矩阵
				for (i=0;i<n;i++){ 
					sigma[i]=Xk1[i]-Xk[i];
					yk[i]=gk1[i]-gk[i];
				}
				
				temp1=0;
				for (i=0;i<n;i++) temp1=temp1+sigma[i]*yk[i];	//计算Mk的分母
				for (i=0;i<n;i++) {
					for (j=0;j<n;j++) {
						Mk[i][j]=sigma[i]*sigma[j]/temp1;		//计算Mk
						temp7[i][j]=0;
						temp8[i][j]=0;
					}
					temp4[i]=0;
					temp5[i]=0;
					temp6[i]=0;
				}


				for (i=0;i<n;i++) 
					for (j=0;j<n;j++) {
						temp4[i]=temp4[i]+yk[j]*Ak[j][i];			//Nk的分母前两项相乘
						temp5[i]=temp5[i]+Ak[i][j]*yk[j];			//Nk的分子前两项相乘
					}
					
				for (i=0;i<n;i++)									//Nk分子的前三项相乘的结果
					for (j=0;j<n;j++) temp7[i][j]=temp5[i]*yk[j];
				
					
				temp1=0;
				for (i=0;i<n;i++)
					temp1=temp1+temp4[i]*yk[i];						//Nk的分母
				
				for (i=0;i<n;i++)
					for (j=0;j<n;j++)								//计算Nk的分子
						for (t=0;t<n;t++) temp8[i][j]=temp8[i][j]+temp7[i][t]*Ak[t][j];

				for (i=0;i<n;i++) {
					for (j=0;j<n;j++) {		//计算Nk和Ak
						Nk[i][j]=temp8[i][j]/temp1;
						Ak[i][j]=Ak[i][j]+Mk[i][j]-Nk[i][j];
					}
					Xk[i]=Xk1[i];			//用新计算得到的Xk1对Xk重新赋值
				}
				calculate_g(Xk1, gk);		//DFP 计算后重新计算梯度
			}
			cout<<endl;
		}
		if (w==-1) break;					//如果满足精度条件，则结束迭代
		cout<<endl;
	}

}