bfgs1.cpp

bfgs算法求全域最小数值点

#include<math.h>
#include<iostream.h>
#include<fstream.h>
#include<iomanip.h>
const int length=10;
const double e=0.00001;
const int iitem2 =100 ; //100
const double cgold = 0.381966;  
const double zeps = 0.0000000001;
const double gold=1.618034;
const double limg=100;
const double tiny=1e-20;
const double tol= 0.0001;
const double iitem=200; 
const double eps=0.0000000001;//

//////////////variable setting////////////////

typedef  double array2[length+1][length+1];
typedef  double array1[length+1];

 int n,m,nnc,kk;
 double f;
 double xishu[length+1][length+1];
 array1 pp,xxi,xx0;

////////////////////read in//////////////////////////
void readin()
{  
     
    double f;
	ifstream fin("in.txt",ios::nocreate);
	fin>>n>>m;
/*in.txt contain n represent number of the variable  
 m represent highest level klof  the variable */
	           
   int i,j;
   for (i=1;i<=n;i++) 
	   fin>>xx0[i]; 
   for (i=1;i<=n;i++)
	  {
          for (j=1;j<=m+1;j++)
		  {  
			  fin>>f; 
		      xishu[i][j]=f;
		  };
   };
          
};

/////////////////求初始的函数值//////////////////
double func2(array1 xx,int n)
 
{  
	int i;int j;double s;
   s =0;
   for (i=1;i<=n;i++)
     for (j=1;j<=m+1;j++)
       s =s+xishu[i][j]*pow(xx[i],n+1-j);
   return s;
 };

double transl(double x)
 {
int j ;
array1 xt;
for (j=1;j<=nnc;j++)
      xt[j] =pp[j] + x * xxi[j];
    return func2(xt,nnc);
};

double func(double x)
{
   return transl(x);
};
/////////////////求初始的导数值////////////////////
void dfunc(array1 xx,array1 df)//procedure dfunc(array1 xx;var df array1);
{
array1 s; int i;int j;double hh;

   for (i=1;i<=n;i++)
   {
      s[i]=0;
      for (j=1;j<=m;j++)
	  {
		 hh=pow(xx[i],m-j);
		  s[i] =s[i]+xishu[i][j]*(m+1-j)*hh;
	  };
      df[i] =s[i];
   };
  };
/////////////////function"sign"////////////////////////
int sign(double x)
{
	int t=1;
 if (x<0) t=-1;
 return t;
}
/////////////////////////////////////////////////
double golden(double ax,double bx,double cx,double tol, double* xmin)
{    
	int teme;int kk;
     double d,fu,r,q,p,xm,r1,r2,a,b;
     double u,tempe,dm,v,w,x,ee,fx,fv1,fw;

    a =ax;
    if (cx<ax) a =cx;
    b  = ax;
    if (cx>ax) b =cx;
    v =bx;
    w =v;
    x =v;
    ee =0.0;
    fx =func(x);
    fv1 =fx;
    fw =fx;
    for (kk=1;kk<=iitem2;kk++)
	{
          xm =0.5*(a+b);
          r1 =tol*fabs(x)+zeps;
          r2 =2.0*r1;
          if (fabs(x-xm)<=r2-0.5*(b-a))  break;
          teme =-1;
          if (fabs(ee)>r1)   
           {
             r =(x-w)*(fx-fv1);
             q =(x-v)*(fx-fw);
             p =(x-v)*q-(x-w)*r;
             q =2.0*(q-r);
             if (q>0.0)    p =-p;
             q = fabs(q);
             tempe =ee;
             ee =d;
             dm = fabs(0.5*q*tempe);
             if (( fabs(p)<dm ) && (p>q*(a-x)) && (p<q*(b-x)))   
               {
                 d =p/q;
                 u =x+d;
                 if ((u-a<r2)||(b-u<r2))
                   d = fabs(r1)*sign(xm - x);
                 teme =0;
               };
           };
        if (teme!=0)   
          {
            if (x>= xm)   
              ee =a-x;
            else
              ee =b-x;
            d =cgold*ee;
          };
        if ( fabs(d)>= r1)   
          u =x+d;
        else
          u =x+ fabs(r1)*sign(d);
        fu =func(u);
        if (fu<=fx)   
          {
            if (u>=x)   
              a =x;
            else
              b =x;
            v =w;
            fv1 =fw;
            w =x;
            fw =fx;
            x =u;
            fx =fu;
          }
        else
          {
            if (u<x)  
             a =u;
            else
             b =u;
            if ((fu<=fw)||(w=x))    
              {
                v =w;
                fv1 =fw;
                w =u;
                fw =fu;
              }
            else
              {
                if ((fu<=fv1)||(v=x)||(v=w))   
                  {
                    v =u;
                    fv1 =fu;
                  };
              };
          };
    };
   *xmin =x;
   return fx;
};
////////////////////////
double ax,bx,cx,fa,fb,fc;
////////////////////////
void searchRange()
 {
double r,q,dm,u,ylm,fu;
    fa =func(ax);
    fb =func(bx);
    if (fb>fa)   
      {
        dm =ax; ax =bx; bx =dm;
        dm =fb; fb =fa; fa =dm;
      };
    cx =bx+gold*(bx-ax);
    fc =func(cx);
    while (fb>=fc) 
	{
          r =(bx-ax)*(fb-fc);
          q =(bx-cx)*(fb-fa);
          dm =q-r;
          if ( fabs(dm)<tiny)  dm =tiny;

          u =bx-((bx-cx)*q-(bx-ax)*r)/(2*dm);
          ylm =bx+limg*(cx-bx);
          if ((bx-u)*(u-cx)>0)   
           {
              fu =func(u);
              if (fu<fc)   
               {
                  ax=bx;  bx=u;
                  fa=fb;  fb=fu;
                  goto EXIT;
               }
              else
               {
                  if (fu>fb)   
                    {
                      cx =u;
                      fc =fu;
                     goto EXIT;
                    };
               };
              u =cx+gold*(cx-bx);
              fu =func(u);
           }
         else
            {
              if ((cx-u)*(u-ylm)>0)   
                {
                  fu =func(u);
                  if (fu<fc)    
                    {
                      bx =cx;
                      cx =u;
                      u =cx+gold*(cx-bx);
                      fb =fc;
                      fc =fu;
                      fu =func(u);
                    };
                 }
                else
                    {
                      if ((u-ylm)*(ylm-cx)>=0)   
                        {
                          u =ylm;
                          fu =func(u);
                        }
                     else
                        {
                          u =cx+gold*(cx-bx);
                          fu =func(u);
                        };
                    };
             };
		ax =bx;
		bx =cx;
		cx =u;
		fa =fb;
		fb =fc;
		fc =fu;
	};
EXIT:
cout<<" ";
};
////////////定方向的线性搜索//////////////////////
void lineSearch(array1 p,array1 xi,int n,double* fret)
{ 
int j;
double xmin,xx,fx;
   
    nnc =n;
    for (j=1;j<=n;j++)
      {
        pp[j] =p[j];
        xxi[j] =xi[j];
      };
        ax =0.0;
		bx=1.0;
        searchRange();
	    *fret =golden(ax,bx,cx,tol,&xmin);
    for (j=1;j<=n;j++)
      {
        xi[j] =xmin*xi[j];
        p[j] = p[j]+xi[j];
      };
};

  void DFPmin(array1 p,int n,double ftol, int *kk1,double* f1)
 {
int i,j,its,kk;
double fp,fac,fad,fae,aaa,bbb,f;
array2 hessin;
array1 xi,g,dg,hdg;
f=*f1; kk=*kk1; 

    fp =func2(p,n);
    dfunc(p,g);
    for (i=1;i<=n;i++)
      {
       for (j=1;j<=n;j++)
         hessin[i][j] =0.0;
       hessin[i][i] =1.0;
       xi[i] =-g[i];
      };
    for (its=1;its<=iitem;its++)
      {
         kk =its;
         lineSearch(p,xi,n,&f);
         if (2.0* fabs(f-fp)<=ftol*( fabs(f)+ fabs(fp)+eps))   
           goto exit2;
         fp =f;
         for (i=1;i<=n;i++)
           dg[i] =g[i];
         f =func2(p,n);
         dfunc(p,g);
         for (i=1;i<=n;i++)
           dg[i] =g[i]-dg[i];
         for (i=1;i<=n;i++)
           {
             hdg[i] =0.0;
             for (j=1;j<=n;j++)
                 hdg[i] =hdg[i]+hessin[i][j]*dg[j];
           };
         fac =0.0;
         fae =0.0;
         for (i=1;i<=n;i++)
           {
             fac =fac+dg[i]*xi[i];
             fae =fae+dg[i]*hdg[i];
           };
         fac =1.0/fac;
         fad =1.0/fae;
         for (i=1;i<=n;i++)
           dg[i] =fac*xi[i]-fad*hdg[i];
         for (i=1;i<=n;i++)
           {
             for (j=1;j<=n;j++)
               {
                 aaa =fac*xi[i]*xi[j]-fad*hdg[i]*hdg[j];
                 bbb =fae*dg[i]*dg[j];
                 hessin[i][j] =hessin[i][j]+aaa;//+bbb;
               };
           };
         for (i=1;i<=n;i++)
           {
             xi[i] =0.0;
             for (j=1;j<=n;j++)
              xi[i] =xi[i]-hessin[i][j]*g[j];
           };
      };
exit2:
	*kk1=kk; 
	*f1=f;

};
  void output(array1 ans,int n)
{
  //ofstream fout("out.txt");
  cout<<endl<<"-----------ans------------"<<endl; 
  for (int i=1;i<=n;i++)
  {   
	  if (fabs(ans[i])<tol) ans[i]=0;
	  cout<<setprecision(3)<<setiosflags(ios::fixed)<<ans[i]<<"  "; 
  }
	  cout<<endl<<"minimor answer is"<<f<<endl;
}
void main()
{
  readin();
  DFPmin(xx0,n,tol,&kk,&f);
  output(xx0,n);
};