Powell法求多元函数极小值(要求目标函数有连续二阶导数)
#include <conio.h>#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define XuYimin_Debug_1
#define N 3 /*N维空间*/
double xq[N]={0,0,0};
double fun(double x[N])
{
double y;
/*目标函数f(X)=
(x1-4)^2+(x2-3)^2+1.5*x3^2+(x1-x2+x3-1)^2+sin(x1+x2+x3-7)*sin(x1+x2+x3-7)
可根据需要修改,但需要同时修改维数N值*/
y=pow((x[0]-4),2)+pow((x[1]-3),2)+1.5*x[2]*x[2]
+pow(x[0]-x[1]+x[2]-1,2)+pow(sin(x[0]+x[1]+x[2]-7),2);
return y;
}
double fd(double x[N],double d[N],double alf)
{
int i;
double y,t[N];
for(i=0;i<N;i++)
t[i]=x[i]+alf*d[i];
y=fun(t);
return y;
}
double alf(double x[N],double d[N])
{
double h0=0.1,a0=0.0,ep1;
double a1,a2,a3,y1,y2,y3,c1,c2,ap,yp,ax,yx,h;
ep1=pow(0.1,12);
a1=a0,h=h0;a2=h,y1=fd(x,d,a1),y2=fd(x,d,a2);
if(y2>y1)
{
h=-1*h,a3=a1,y3=y1;
loop1:a1=a2,y1=y2,a2=a3,y2=y3;
loop2:a3=a2+h,y3=fd(x,d,a3); if(y3<y2) {h=2*h;goto loop1;}
}
else goto loop2;
if( fabs(y2)<ep1 ) goto loop5;
loop3:
c1=(y3-y1)/(a3-a1);
c2=((y2-y1)/(a2-a1)-c1)/(a2-a3);
if(c2==0) goto loop4;
ap=0.5*(a1+a3-c1/c2);
yp=fd(x,d,ap);
while(!fabs((y2-yp)/y2)<ep1)
{
if((ap-a2)*h>0)
{
if(y2>=yp) a1=a2,y1=y2,a2=ap,y2=yp;
else a3=ap,y3=yp;
}
else
{
if(y2>=yp) a3=a2,y3=y2,a2=ap,y2=yp;
else a1=ap,y1=yp;
}
goto loop3;
}
if(y2<yp) ax=a2,yx=y2;
else ax=ap,yx=yp;
return ax;
loop4:return a2;
loop5:return 0;
}
double maxt(double del[N+1])
{
int i;
double y=del[0];
for (i=1;i<N+1;i++)
if(del[i]>y)
y=del[i];
return y;
}
void powell()
{
int i,j,id,pd[3],xym=1,m;
double al,esp,ALF[N+1],x[2],pa[N+1],X[N][N+1],del[N+1],d[N][N+1],
F[3],fu[N+1],te[N],tem[N],dk[N],xk[N],delm;
esp=pow(0.1,10);x[0]=xq[0];x[1]=xq[1];
for(i=0;i<N;i++)
for(j=0;j<N+1;j++)
d[i][j]=0;
for(i=0;i<N;i++) d[i][i+1]=1;
loop:
#ifdef XuYimin_Debug_1
printf("fangxiang juzhen:d=\n");
for(i=0;i<N;i++)
for(j=0;j<N+1;j++)
{
printf("%15g",d[i][j]);
if(j==N) printf("\n");
}
printf("\n");
#endif
for(i=0;i<N;i++) X[i][0]=x[i];
ALF[0]=0;del[0]=0;
fu[0]=fun(x);
for(id=1;id<N+1;id++)
{
for(i=0;i<N;i++) tem[i]=d[i][id];
ALF[id]=alf(x,tem);
for(i=0;i<N;i++) te[i]=x[i]+ALF[id]*d[i][id];
fu[id]=fun(te);
for(i=0;i<N;i++) X[i][id]=te[i];
for(i=0;i<N;i++) x[i]=te[i];
}
for(id=1;id<N+1;id++) del[id]=fu[id-1]-fu[id];
#ifdef XuYimin_Debug_1
printf("zuijia buchang juzhen:ALF:\n");
for(i=0;i<N+1;i++) printf("%15g",ALF[i]);
printf("\n\ndiedaidian juzhen:X=\n");
for(i=0;i<N;i++)
for(id=0;id<N+1;id++)
{
printf("%15g",X[i][id]);
if(id==N) printf("\n");
}
printf("\nhanshuzhi juzhen:fu=\n");
for(id=0;id<N+1;id++) printf("%15g",fu[id]);
printf("\n\nhanshuzhi chazhi juzhen:del=\n");
for(id=0;id<N+1;id++) printf("%15g",del[id]);
printf("\n");
#endif
for(i=0;i<N;i++)
{
dk[i]=X[i][N]-X[i][0];
xk[i]=2*X[i][N]-X[i][0];
}
for(i=0;i<N;i++) x[i]=X[i][0];
F[0]=fun(x);
for(i=0;i<N;i++) x[i]=X[i][N];
F[1]=fun(x);
F[2]=fun(xk);
delm=maxt(del);
for(i=0;i<N+1;i++)
if(del[i]==delm)
m=i;
#ifdef XuYimin_Debug_1
printf("x0_k=\n");
for(i=0;i<N;i++) printf("%15g\n",X[i][0]);
printf("\nxn_k=\n");
for(i=0;i<N;i++) printf("%15g\n",X[i][N]);
printf("\nxn+1_k=\n");
for(i=0;i<N;i++) printf("%15g\n",xk[i]);
printf("\nhanshuzhi zengliang juzhen F=\n");
for(i=0;i<3;i++) printf("%15g",F[i]);
printf("\n\nzuida chazhi:m=%d delm=%g\n\n",m,delm);
#endif
pd[0]=(F[2]<F[0]);
pd[1]=((F[0]-2*F[1]+F[2])*(F[0]-F[1]-delm)*(F[0]-F[1]-delm)
<0.5*delm*(F[0]-F[2])*(F[0]-F[2]));
pd[2]=(F[1]<F[2]);
#ifdef XuYimin_Debug_1
printf("pd[0]=%d pd[1]=%d pd[2]=%d\n\n",pd[0],pd[1],pd[2]);
#endif
if(pd[0]*pd[1]==1)
{
for(i=0;i<N;i++) te[i]=X[i][N];
al=alf(te,dk);
for(i=0;i<N;i++) x[i]=X[i][N]+al*dk[i];
for(j=m;j<N;j++)
for(i=0;i<N;i++)
d[i][j]=d[i][j+1];
for(i=0;i<N;i++) d[i][N]=dk[i];
#ifdef XuYimin_Debug_1
printf("xin diedaidian:x=\n");
for(i=0;i<N;i++) printf("%15g\n",x[i]);
printf("\nxin fangxiang juzhen:d=\n");
for(i=0;i<N;i++)
for(j=0;j<N+1;j++)
{
printf("%15g",d[i][j]);
if(j==N) printf("\n");
}
printf("\n");
#endif
}
else
if(pd[2]==1)
for(i=0;i<N;i++) x[i]=X[i][N];
else
for(i=0;i<N;i++) x[i]=xk[i];
pa[N]=0;
for(i=0;i<N;i++) pa[i]=X[i][N]-X[i][0];
for(i=0;i<N;i++) pa[N]+=pa[i]*pa[i];
pa[N]=sqrt(pa[N]);
#ifdef XuYimin_Debug_1
printf("piancha juzhen:PA=\n");
for(i=0;i<N;i++)
printf("%15g\n",pa[i]);
printf("\npiancha juli:pa=%g\n\n",pa[N]);
printf("************Di %d ci diedai jieshu!************\n\n\n\n\n",xym++);
#endif
if(pa[N]<esp)
{
printf("The Result:\n");
printf("jizhi dian:X=\n");
for(i=0;i<N;i++)
printf("%15g\n",x[i]);
printf("\njizhidian:f(X)=%f",fun(x));
}
else goto loop;
}
int main()
{
powell();
getch();
return 0;
}
程序刚编完,正在修改,希望大家提点意见!