hbfgsopt.h

《VC++ 编程技巧与示例 .rar》各个示例代码绝对可用

//类模板使用示例：
// HDFPOpt.h: interface for the CBFGSOpt class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_HDFPOPT_H__18F1F1DA_10A2_11D0_9458_0000000038B2__INCLUDED_)
#define AFX_HDFPOPT_H__18F1F1DA_10A2_11D0_9458_0000000038B2__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "CMatrix.h"
#include "math.h"
#define ALFA  1.61803398875
#define GOLD  0.61803398875
#define MAXI  150//最多迭代次数；
template <class T>
class CBFGSOpt  
{
public:
	double * X;//设计变量；
	double F[MAXI];//目标函数值；
	double fhs[500];//记录迭代历史的数组；
	int      N;//设计变量数目；
	BOOL Inited;//是否已经初始化了；
	BOOL Success;//优化是否成功；
	T * t;
	double (T::* f    )(double x0[50],double y0[50]);//目标函数；

	//与寻优有关的变量：
	CMatrix A ;//A矩阵；
	CMatrix S ;//搜索方向；
	CMatrix df;//f的梯度矩阵；
	CMatrix df1;//f的梯度矩阵；
	CMatrix dltx;
	CMatrix dltg;
	double theta;
	double Lmd;//搜索步长；
	int count; //迭代次数；

public:
	CBFGSOpt()
	{
		Inited=FALSE;
		X=NULL;
		f=NULL;
		Success=FALSE;
	}

	virtual ~CBFGSOpt()
	{
		if(X!=NULL)
		{
			delete [] X;
			X=NULL;
		}
	}

	void InitData(T *tt,int n)
	{
		t=tt;
		N=n;
		if(X==NULL)
		{
			X=new double [N];
		}
		else
		{
			delete [] X;
			X=new double [N];
		}

		//初始化矩阵：
		A   .setvalue(N,N);//A矩阵；
		S   .setvalue(N,1);//搜索方向；
		df  .setvalue(N,1);//f的梯度矩阵；
		df1 .setvalue(N,1);//f的梯度矩阵；
		dltx.setvalue(N,1);
		dltg.setvalue(N,1);

		if(f!=0)
		{
			Inited=TRUE;
		}
	}

	double CalF(double x[100],double y[100])
	{
		double re;
		re=(t->*f)(x,y);
		return re;
	}

	double FindMin(double a ,double b ,double c ,
				   double fa,double fb,double fc)
	{
		double re;
		double fz,fm;

		fz=(b*b-c*c)*fa+(c*c-a*a)*fb+(a*a-b*b)*fc;
		fm=(b-c)*fa+(c-a)*fb+(a-b)*fc;

		re=1/2*fz/fm;
		return re;
	}

	double CalFLmd(double lmd,CMatrix s,
				   double x[100],double y[100])
	{
		double re;
		int i;
		for(i=0;i<N;i++)
		{
			x[i]=x[i]+lmd*s(i,0);
		}
		re=CalF(x,y);
		for(i=0;i<N;i++)
		{
			x[i]=x[i]-lmd*s(i,0);
		}
		return re;
	}

	void FindStepRange(
						double &ll   ,double &lu,
						double x[100],double y[100]
					   )
	{
		double fl,fu,fz;
		double lz;
		double lmax;
		lmax=1e20;
		ll=0.0;
		lu=1.0;
		fl=CalFLmd(ll,S,x,y);
		fu=CalFLmd(lu,S,x,y);
		if(fu>fl)
		{
			return;
		}
		lz=lu;
		fz=fu;
		lu=(1+ALFA)*lz-ALFA*ll;
		
		if(lu>lmax)
		{
			AfxMessageBox("The minimum is un-bouned!");
			count=MAXI;
			return;
		}
		fu=CalFLmd(lu,S,x,y);
		if(lu>lz)
		{
			return;
		}
		ll=lz;
		fl=fz;
	}

	double FindStep(double x[100],double y[100])
	{
		double re;
/*
	问题为：
	min f(Lmd)
*/	
		double al,b,c,au;
		double fal,fb,fc,fau;
		double d;
		double dlt;
		FindStepRange(al,au,x,y);
		d=au-al;
		if(d<0)
		{
			AfxMessageBox("d<0!");
			count=MAXI;
			return re;
		}
		b=al+(1-GOLD)*d;
		c=al+GOLD*d;
		dlt=1000.0;
		while(dlt>1e-8)
		{
	
			fal=CalFLmd(al,S,x,y);
			fb =CalFLmd(b ,S,x,y);
			fc =CalFLmd(c ,S,x,y);
			fau=CalFLmd(au,S,x,y);

			if(fb>=fc)//b will be on the bound;
			{
				al=b;
				b=c;
				d=au-al;
				c=al+GOLD*d;
			}
			else//c will be on the bound;
			{
				au=c;
				c=b;
				d=au-al;
				b=al+(1-GOLD)*d;
			}
			dlt=fabs(c-b);
		}
		re=c;
		return re;
	}

 CMatrix Finddf(double x[100],double y[100])
	{
		int i;
		CMatrix ddf;
		double * df0;
		double f0,f1,f2;
		double step,h;
		df0=new double[N];
		step=1.0e5;
		f0=CalF(x,y);
		//按设计变量循环：
		for(i=0;i<N;i++)
		{
			if(x[i]!=0.0)//x[i]不为0.0;
			{
				h=x[i]/step;
			}
			else
			{
				h=1/step;
			}
			x[i]=x[i]+h;
			f1=CalF(x,y);
			x[i]=x[i]+h;
			f2=CalF(x,y);
			x[i]=x[i]-2*h;
			df0[i]=(-3*f0+4*f1-f2)/2.0/h;
		}
		ddf.setvalue (N,1,df0);
		delete [] df0;
		fhs[count]=f1;
		return ddf;
	}

	double FindTheOpt(double x[100],double y[100])
	{
		double re;
		double modf;
		int i;
		count=0;
		if(Inited==FALSE)
		{
			re=0;
			return re;
		}
		//如果count==0，A=I：
		if(count==0)
		{
			for(i=0;i<N;i++)
			{
				A.setelem(i,i,1.0);
			}
		}
		modf=1000.0;
		Success=FALSE;
		while((count<MAXI)&&(Success==FALSE))
		{
			df=Finddf(x,y);
			S=-A*df;
			Lmd=FindStep(x,y);
			for(i=0;i<N;i++)
			{
				x[i]=x[i]+S(i,0)*Lmd;
			}
			df1=Finddf(x,y);
			re=CalF(x,y);
			F[count]=re;
			modf=df1.mod();
			if(modf<0.0001)
			{
				Success=TRUE;
				continue;
			}
			//求新的A：
			dltx=Lmd*S;
			dltg=df1-df;
			double  temp1,temp2,temp3;
			CMatrix t1(1,1);
			//分子：
			t1=(~dltx)*dltg;
			temp1=t1(0,0);
			//θ分子：
			t1=(~dltg)*A*dltg;
			temp2=t1(0,0);
			//θ分母：
			temp3=temp1;
			//θ：
			theta=1+temp2/temp3;
			A=A+1/temp1*(
						theta*dltx*(~dltx)-
						A*dltg*(~dltx)-
						dltx*(~dltg)*A
						);
			count++;
		}
		//记录最优解：
		for(i=0;i<N;i++)
		{
			X[i]=x[i];
		}
		return re;
	}

	BOOL SaveInfomation(CString fn)
	{
		int i;
		FILE *fp;
		if((fp=fopen(fn,"wt"))==NULL)
		{
			AfxMessageBox("BFGSOpt 类存盘错误！");
			return FALSE;
		}
		fprintf(fp,"\n\n\t☆☆☆☆☆☆☆☆☆☆☆☆☆\n\n");
		fprintf(fp,"\tBFGSOpt 优化结果：\n\n");
		fprintf(fp,"\n\n\t☆☆☆☆☆☆☆☆☆☆☆☆☆\n\n");
		
		fprintf(fp,"\n\n\t迭代次数：\n\n\t");
		fprintf(fp,"\n\n\tcount=%d\n\n\t",count);

		fprintf(fp,"\n\n\t目标函数值：\n\n\t");
		for(i=0;i<=count;i++)
		{
			fprintf(fp,"%10.4g     ",F[i]);
			if(i%5==0)
			{
				fprintf(fp,"\n\t");
			}
		}
		fprintf(fp,"\n\n\t设计变量最优解：\n\n\t");
		for(i=0;i<N;i++)
		{
			fprintf(fp,"%10.4g    ",X[i]);
			if(i%5==0)
			{
				fprintf(fp,"\n\t");
			}
		}
		fclose(fp);	
		fprintf(fp,"\n\n\t☆☆☆☆☆☆☆☆☆☆☆☆☆\n\n");
	}
};

#endif // !defined(AFX_HDFPOPT_H__18F1F1DA_10A2_11D0_9458_0000000038B2__INCLUDED_)