机械优化设计大作业：平面连杆机构的优化设计，用C语言编程！

计算 f(x1,x2)=x1^2+2*x2^2-4*x1-2*x1*x2 的无约束极值，初始点x0=[1,1]。

/*
tt ---- 一维搜索初始步长
ff ---- 差分法求梯度时的步长
ac ---- 终止迭代收敛精度
ad ---- 一维搜索收敛精度
n ----- 设计变量的维数
xk[n] -- 迭代初始点
*/

#include
#include
#include
#include

#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 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] 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 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 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{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{
s[i]=0;
for(j=0;j {
h[i][j]=0;
if(j==i) h[i][j]=1;
}
}
g1=gradient(xk);
fi=fny(xk);
x0=xk;

for(k=0;k {
ib=0;
if(ia==1)
ib=0;
for(i=0;i for(i=0;i for(j=0;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 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 for(j=0;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 {
dg[i]=g2[i]-g1[i];
dx[i]=xk1[i]-xk[i];
}

for(i=0;i {
int j;
u[i]=0;v[i]=0;
for(j=0;j {
u[i]=u[i]+dg[j]*h[j][i];
v[i]=v[i]+dg[j]*h[i][j];
}
}

for(j=0;j {
a1+=dx[j]*dg[j];
a2+=v[j]*dg[j];
}
if(fabs(a1)!=0)
{
a2=1+a2/a1;
for(i=0;i for(j=0;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
printf("\n\tf*=%f",f);
getch();
}

这是基于一本书上的算法。但我很奇怪，原书中的算法有结果列出，但是我却不能通过编译。真是纳闷！修改后可以得到结果了，如果你要使用这个简单的程序，你只需更改 维数n、double fny(double *x)的实现部分以及main函数中的xk初值就可以了。不过这个程序也不是很好。