powell.cpp

鲍威尔法无约束优化程序

/*本程序包含5个C文件: mpowell.c,  powell.c,
funct.c(目标函数), jtf.c(进退法), hjfgf.c(黄金分割法)*/
//题目 y=x1*x1+x2*x2-x1*x2-10*x1-4*x2+60

//#include "powell.c"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#define pi  3.1415926/4

double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[]);
double powell(double p[],double h0,double eps,double epsg,int n,double x[]);
double objf(double x[]);
void jtf(double x0[],double h0,double s[],int n,double a[],double b[]);
double gold(double a[],double b[],double eps,int n,double xx[]);
void main()
{double p[]={23,0.8,42,12,16};
double ff,x[5];
ff=powell(p,0.3,0.0001,0.00001,5,x);
printf("f0=%F\n",objf(p));
printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff);
 
}
//powell.c代码如下： 
//#include "hjfgf.c"
//一维优化求极小点函数oneoptim
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
{double *a,*b,ff;
a=(double *)malloc(n*sizeof(double));
b=(double *)malloc(n*sizeof(double));
jtf(x0,h0,s,n,a,b);
ff=gold(a,b,epsg,n,x);
free(a);
free(b);
return (ff);
}
//powell法程序
double powell(double p[],double h0,double eps,double epsg,int n,double x[])
{int i,j,m;
double *xx[4],*ss,*s;
double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
ss=(double *)malloc(n*(n+1)*sizeof(double));
s=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
  {for(j=0;j<=n;j++)
   *(ss+i*(n+1)+j)=0;
   *(ss+i*(n+1)+i)=1;
  }
  for(i=0;i<4;i++)
  xx[i]=(double *)malloc(n*sizeof(double));
  for(i=0;i<n;i++)
  *(xx[0]+i)=p[i];
  for(;;)
  {for(i=0;i<n;i++)
    {*(xx[1]+i)=*(xx[0]+i);
     x[i]=*(xx[1]+i);
    }
   f0=f1=objf(x);
   dlt=-1;
   for(j=0;j<n;j++)
    {for(i=0;i<n;i++)
      {*(xx[0]+i)=x[i];
       *(s+i)=*(ss+i*(n+1)+j);
      }
     f=oneoptim(xx[0],s,h0,epsg,n,x);
     df=f0-f;
     if(df>dlt)
      {dlt=df;
       m=j;
      }
    }
    sdx=0;
    for(i=0;i<n;i++)
     sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
    if(sdx<eps)
     {free(ss);
      free(s);
      for(i=0;i<4;i++)
      free(xx[i]);
      return(f);
     }
    for(i=0;i<n;i++)
     *(xx[2]+i)=x[i];
    f2=f;
    for(i=0;i<n;i++)
    {*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
     x[i]=*(xx[3]+i);
    }
    fx=objf(x);
    f3=fx;
    q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
    d=0.5*dlt*(f1-f3)*(f1-f3);
    if((f3<f1)||(q<d))
     {if(f2<=f3)
       for(i=0;i<n;i++)
        *(xx[0]+i)=*(xx[2]+i);
      else
       for(i=0;i<n;i++)
        *(xx[0]+i)=*(xx[3]+i);
     }
    else
    {for(i=0;i<n;i++)
     {*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
      *(s+i)=*(ss+(i+1)*(n+1));
     }
   f=oneoptim(xx[0],s,h0,epsg,n,x);
   for(i=0;i<n;i++)
    *(xx[0]+i)=x[i];
   for(j=m+1;j<=n;j++)
   for(i=0;i<n;i++)
    *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
    }
  }
}
//funct.c代码如下： 

double objf(double x[])    //目标函数
{double ff;
//ff=x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60;
 double u1,u2,u3,u4;
u1=x[2]*(x[3]*x[3]+x[4]*x[4])+28*pow(x[3],2)+32*pow(x[4],2);
u2=0.92*x[0]*x[4]*x[4]-x[0]*x[3]*x[3]-85*x[0]*x[1]*x[1];
u3=-(1.6*x[0]*x[1]*x[4]+776.16*x[0]);
u4=1428*x[0]*x[1]+13.44*x[0]*x[4];
ff=pi*(u1+u2+u3+u4);
return(ff);
}
//jtf.c代码如下： 
//#include "funct.c"
void jtf(double x0[],double h0,double s[],int n,double a[],double b[])
{int i;
double *x[3],h,f1,f2,f3;
for(i=0;i<3;i++)
  x[i]=(double *)malloc(n*sizeof(double));
  h=h0;
for(i=0;i<n;i++)
  *(x[0]+i)=x0[i];
f1=objf(x[0]);
for(i=0;i<n;i++)
  *(x[1]+i)=*(x[0]+i)+h*s[i];
  f2=objf(x[1]);
if(f2>=f1)
  {h=-h0;
    for(i=0;i<n;i++)
    *(x[2]+i)=*(x[0]+i);
   f3=f1;
    for(i=0;i<n;i++)
    {*(x[0]+i)=*(x[1]+i);
     *(x[1]+i)=*(x[2]+i);
    }
   f1=f2;
   f2=f3;
   }
   for(;;)
   {h=2*h;
     for(i=0;i<n;i++)
     *(x[2]+i)=*(x[1]+i)+h*s[i];
   f3=objf(x[2]);
   if(f2<f3) break;
   else
    { for(i=0;i<n;i++)
       {*(x[0]+i)=*(x[1]+i);
        *(x[1]+i)=*(x[2]+i);
       }
      f1=f2;
      f2=f3;
    }
   }
   if(h<0)
    for(i=0;i<n;i++)
    {a[i]=*(x[2]+i);
     b[i]=*(x[0]+i);
    }
   else
    for(i=0;i<n;i++)
    {a[i]=*(x[0]+i);
     b[i]=*(x[2]+i);
     }
   for(i=0;i<3;i++)
   free(x[i]);
}
//黄金分割法 hjfgf.c代码如下： 
//#include "jtf.c"
double gold(double a[],double b[],double eps,int n,double xx[])
{int i;
double f1,f2,*x[2],ff,q,w;
for(i=0;i<2;i++)
  x[i]=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
  {*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
   *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
  }
  f1=objf(x[0]);
  f2=objf(x[1]);
  do
   {if(f1>f2)
     {for(i=0;i<n;i++)
      {b[i]=*(x[0]+i);
       *(x[0]+i)=*(x[1]+i);
       }
     f1=f2;
     for(i=0;i<n;i++)
      *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
     f2=objf(x[1]);
     }
    else
     { for(i=0;i<n;i++)
       {a[i]=*(x[1]+i);
       *(x[1]+i)=*(x[0]+i);}
     f2=f1;
    for(i=0;i<n;i++)
     *(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
    f1=objf(x[0]);
     }
  q=0;
  for(i=0;i<n;i++)
   q=q+(b[i]-a[i])*(b[i]-a[i]);
  w=sqrt(q);
  }while(w>eps);
  for(i=0;i<n;i++)
   xx[i]=0.5*(a[i]+b[i]);
  ff=objf(xx);
  for(i=0;i<2;i++)
  free(x[i]);
  return(ff);
}