/*
tt ---- 壹維搜索初始步長
ff ---- 差分法求梯度時的步長
ac ---- 終止叠代收斂精度
ad ---- 壹維搜索收斂精度
n ----- 設計變量的維數
xk[n] -- 叠代初始點
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<conio.h>
#define tt 0.01
#define ff 1.0e-6
#define ac 1.0e-6
#define ad 1.0e-6
#define n 2
double ia;
double fny(double *x)
{
double x1=x[0],x2=x[1];
double f;
f=x1*x1+2*x2*x2-4*x1-2*x1*x2;
return f;
}
double * iterate(double *x,double a,double *s)
{
double *x1;
int i;
x1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
x1[i]=x[i]+a*s[i];
return x1;
}
double func(double *x,double a,double *s)
{
double *x1;
double f;
x1=iterate(x,a,s);
f=fny(x1);
return f;
}
void finding(double a[3],double f[3],double *xk,double *s)
{
double t=tt;
int i;
double a1,f1;
a[0]=0;f[0]=func(xk,a[0],s);
for(i=0;;i++)
{
a[1]=a[0]+t;
f[1]=func(xk,a[1],s);
if(f[1]<f[0]) break;
if(fabs(f[1]-f[0])>=ad)
{
t=-t;
a[0]=a[1];f[0]=f[1];
}
else
{
if(ia==1) return; //break
t=t/2;ia=1;
}
}
for(i=0;;i++)
{
a[2]=a[1]+t;
f[2]=func(xk,a[2],s);
if(f[2]>f[1]) break;
t=2*t;
a[0]=a[1];f[0]=f[1];
a[1]=a[2];f[1]=f[2];
}
if(a[0]>a[2])
{
a1=a[0];
f1=f[0];
a[0]=a[2];
f[0]=f[2];
a[2]=a1;
f[2]=f1;
}
return;
}
double lagrange(double *xk,double *ft,double *s)
{
int i;
double a[3],f[3];
double b,c,d,aa;
finding(a,f,xk,s);
for(i=0;;i++)
{
if(ia==1)
d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2]);
if(fabs(d)==0) break;
c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d;
if(fabs(c)==0) break;
b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1]);
aa=-b/(2*c);
*ft=func(xk,aa,s);
if(fabs(aa-a[1])<=ad)
if(aa>a[1])
{
if(*ft>f[1])
else if(*ft<f[1])
else if(*ft==f[1])
{
a[2]=aa;a[0]=a[1];
f[2]=*ft;f[0]=f[1];
a[1]=(a[0]+a[2])/2;
f[1]=func(xk,a[1],s);
}
}
else
{
if(*ft>f[1])
else if(*ft<f[1])
else if(*ft==f[1])
{a[0]=aa;a[2]=a[1];
f[0]=*ft;f[2]=f[1];
a[1]=(a[0]+a[2])/2;
f[1]=func(xk,a[1],s);
}
}
}
if(*ft>f[1])
return aa;
}
double *gradient(double *xk)
{
double *g,f1,f2,q;
int i;
g=(double*)malloc(n*sizeof(double));
f1=fny(xk);
for(i=0;i<n;i++)
{q=ff;
xk[i]=xk[i]+q; f2=fny(xk);
g[i]=(f2-f1)/q; xk[i]=xk[i]-q;
}
return g;
}
double * bfgs(double *xk)
{
double u[n],v[n],h[n][n],dx[n],dg[n],s[n];
double aa,ib;
double *ft,*xk1,*g1,*g2,*xx,*x0=xk;
double fi;
int i,j,k;
ft=(double *)malloc(sizeof(double));
xk1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{
s[i]=0;
for(j=0;j<n;j++)
{
h[i][j]=0;
if(j==i) h[i][j]=1;
}
}
g1=gradient(xk);
fi=fny(xk);
x0=xk;
for(k=0;k<n;k++)
{
ib=0;
if(ia==1)
ib=0;
for(i=0;i<n;i++) s[i]=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
s[i]+= -h[i][j]*g1[j];
aa=lagrange(xk,ft,s);
xk1=iterate(xk,aa,s);
g2=gradient(xk1);
for(i=0;i<n;i++)
if((fabs(g2[i])>=ac)&&(fabs(g2[i]-g1[i])>=ac))
if(ib==0)
fi=*ft;
if(k==n-1)
{ int j;
xk=xk1;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
h[i][j]=0;
if(j==i) h[i][j]=1;
}
g1=g2; k=-1;
}
else
{
int j;
double a1=0,a2=0;
for(i=0;i<n;i++)
{
dg[i]=g2[i]-g1[i];
dx[i]=xk1[i]-xk[i];
}
for(i=0;i<n;i++)
{
int j;
u[i]=0;v[i]=0;
for(j=0;j<n;j++)
{
u[i]=u[i]+dg[j]*h[j][i];
v[i]=v[i]+dg[j]*h[i][j];
}
}
for(j=0;j<n;j++)
{
a1+=dx[j]*dg[j];
a2+=v[j]*dg[j];
}
if(fabs(a1)!=0)
{
a2=1+a2/a1;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
h[i][j]+=(a2*dx[i]*dx[j]-v[i]*dx[j]-dx[i]*u[j])/a1;
}
xk=xk1; g1=g2;
}
}
if(*ft>fi)
xk=x0;
return xx;
}
void main ()
{
int k;
double *xx,f;
double xk[n]=;
xx=bfgs(xk);
f=fny(xx);
printf("\n\nThe Optimal Design Result Is:\n");
for(k=0;k<n;k++)
printf("\n\tf*=%f",f);
getch();
}
這是基於壹本書上的算法。但我很奇怪，原書中的算法有結果列出，但是我卻不能通過編譯。真是納悶！修改後可以得到結果了，如果妳要使用這個簡單的程序，妳只需更改 維數n、double fny(double *x)的實現部分以及main函數中的xk初值就可以了。不過這個程序也不是很好。