/*
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=0i<ni++)
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]=0f[0]=func(xk,a[0],s)
for(i=0i++)
{
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/2ia=1
}
}
for(i=0i++)
{
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=0i++)
{
if(ia==1) { aa=a[1]*ft=f[1]break}
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(*ft>f[1]) aa=a[1]break}
if(aa>a[1])
{
if(*ft>f[1]) {a[2]=aaf[2]=*ft}
else if(*ft<f[1]) {a[0]=a[1]a[1]=aaf[0]=f[1]f[1]=*ft}
else if(*ft==f[1])
{
a[2]=aaa[0]=a[1]
f[2]=*ftf[0]=f[1]
a[1]=(a[0]+a[2])/2
f[1]=func(xk,a[1],s)
}
}
else
{
if(*ft>f[1]) {a[0]=aaf[0]=*ft}
else if(*ft<f[1]) {a[2]=a[1]a[1]=aaf[2]=f[1]f[1]=*ft}
else if(*ft==f[1])
{a[0]=aaa[2]=a[1]
f[0]=*ftf[2]=f[1]
a[1]=(a[0]+a[2])/2
f[1]=func(xk,a[1],s)
}
}
}
if(*ft>f[1]) {*ft=f[1]aa=a[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=0i<ni++)
{q=ff
xk[i]=xk[i]+qf2=fny(xk)
g[i]=(f2-f1)/qxk[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=0i<ni++)
{
s[i]=0
for(j=0j<nj++)
{
h[i][j]=0
if(j==i) h[i][j]=1
}
}
g1=gradient(xk)
fi=fny(xk)
x0=xk
for(k=0k<nk++)
{
ib=0
if(ia==1) { xx=xkbreak}
ib=0
for(i=0i<ni++) s[i]=0
for(i=0i<ni++)
for(j=0j<nj++)
s[i]+= -h[i][j]*g1[j]
aa=lagrange(xk,ft,s)
xk1=iterate(xk,aa,s)
g2=gradient(xk1)
for(i=0i<ni++)
if((fabs(g2[i])>=ac)&&(fabs(g2[i]-g1[i])>=ac))
{ib=ib+1}
if(ib==0) { xx=xk1break}
fi=*ft
if(k==n-1)
{ int j
xk=xk1
for(i=0i<ni++)
for(j=0j<nj++)
{
h[i][j]=0
if(j==i) h[i][j]=1
}
g1=g2k=-1
}
else
{
int j
double a1=0,a2=0
for(i=0i<ni++)
{
dg[i]=g2[i]-g1[i]
dx[i]=xk1[i]-xk[i]
}
for(i=0i<ni++)
{
int j
u[i]=0v[i]=0
for(j=0j<nj++)
{
u[i]=u[i]+dg[j]*h[j][i]
v[i]=v[i]+dg[j]*h[i][j]
}
}
for(j=0j<nj++)
{
a1+=dx[j]*dg[j]
a2+=v[j]*dg[j]
}
if(fabs(a1)!=0)
{
a2=1+a2/a1
for(i=0i<ni++)
for(j=0j<nj++)
h[i][j]+=(a2*dx[i]*dx[j]-v[i]*dx[j]-dx[i]*u[j])/a1
}
xk=xk1g1=g2
}
}
if(*ft>fi) { *ft=fixx=xk}
xk=x0
return xx
}
void main ()
{
int k
double *xx,f
double xk[n]={1,1}
xx=bfgs(xk)
f=fny(xx)
printf("\n\nThe Optimal Design Result Is:\n")
for(k=0k<nk++)
{printf("\n\tx[%d]*=%f",k+1,xx[k])}
printf("\n\tf*=%f",f)
getch()
}
这是基于一本书上的算法。但我很奇怪，原书中的算法有结果列出，但是我却不能通过编译。真是纳闷！修改后可以得到结果了，如果你要使用这个简单的程序，你只需更改 维数n、double fny(double *x)的实现部分以及main函数中的xk初值就可以了。不过这个程序也不是很好。
#include <stdio.h>int f(int a,int b,int c,int d){
return abs((a+b)*(c+d)+2)
}
int main()
{
int a,b,c,d
printf("please enter a,b,c,d\n")
scanf("%d,%d,%d,%d",&a,&b,&c,&d)
printf("result=%d\n",f(a,b,c,d))
getch()
}
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）