mpowell.c.txt

优化算法,POWELL算法 C码,不需要目标函数求导,适应于多变量优化.

mpowell.c代码如下：

#include "powell.c"
main()
{double p[]={1,2};
double ff,x[2];
ff=powell(p,0.3,0.001,0.0001,2,x);
  printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff);
  getch();
}
powell.c代码如下：
#include "hjfgf.c"
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);
}
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代码如下：
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
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;
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);
}