powell1.c

powell优化算法的C语言原代码

/*powell算法程序*/
/*目标函数*/
#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);
 }

 /*进退法子程序，用于确定函数初始搜索区间*/
 /*输出变量：a:n维数组，搜索区间的下限*/
 /*          b:n维数组，搜索区间的下限*/
 /*输入变量：xo[]:n维数组的搜索初始点*/
 /*          h0:初始步长            */
 /*          s:n维数组的搜索初始点  */
 /*          n:优化模型的维数       */

 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));
 getch();
  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);
     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<n;i++)     /*释放内存空间，完成搜索*/
     free(x[i]);
     }

     /*黄金分割法子程序用来解决一维优化问题   */
     /*输出变量：xx:最优点，程序返回最优目标值*/
     /*输入变量: a: 寻优区间的下限            */
     /*          b: 寻优区间的上限            */
     /*          eps:迭代精度                 */
     /*          n:优化模型的维数             */
 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<n;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);
     }

   /*用来解N维无约束优化问题*/
 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);
  }
/*输出变量:x:目标函数最优点，程序返回最优目标值*/
/*输入变量:p:寻优初始点*/
/*         h0:进退法寻优迭代精度*/
/*         eps:总体寻优迭代精度*/
/*         N:优化模型维数*/
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));
 printf("step into powell\n");
 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];
                                                /*把初始点赋值给X变量*/
  for(;;)                                      /*无限循环，通过return返回*/
 {
   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];                 /*将寻优获得的迭代值赋值给XX变量*/
     *(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);                   /*跳出循环体，POWELL结束*/
        }
     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); /*powell条件*/
     d=0.5*dlt*(f1-f3)*(f1-f3);
     /*如果至少有一个条件满足时，搜索方向不变，按原来的所搜方向。在此条件下
       如果f2>f3则最好一次最优方向作为下一轮迭代的初始点，否则，附件方向
      xx[3]+i =2.*(*xx[2]+i)-(*(xx[1]+i)) 作为初始点*/
     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);
           }         /*else对应的结束符*/
         }           /*无限循环结束符*/
    }                /*函数结束符*/
 /*主程序，求函数局部最优值，输出最优点及其目标函数值*/
 void main()
 {double p[]={1,2};
  double ff,x[2];
  ff=powell(p,0.3,0.001,0.0001,2,x);
  printf("put out the value of object\n");
  printf("x[0]=%f,x[1]=%f,ff=%f",x[0],x[1],ff );
  getch();
  }