够难！帮忙处理一下floating point error :domain(急用)
出错处已经标出#include "Stdio.h"
#include "Conio.h"
#include "stdlib.h"
#define pi 3.1415
#define kkg 9
double r0=1;
double Dealpos1(int i)
{double a[5]={ 14.5,15.8,16,15,12.5};
return a[i];
}
double Dealpos2(int i)
{double b[5]={1.8,6,9.5,14.5,19};
return b[i];
}
double objf(double x[])
{float f=0;
float Lmin=1.772,dlt_L=0.25;
float c1,c2,c3,L,d1,d2,aa,bb;
int n;
for(n=0;n<5;n++)
{ L=Lmin+dlt_L*n;
c1=acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]));
c2=acos(x[5]-x[0]*cos(x[1]))/sqrt(x[0]*x[0]+x[5]*x[5]-2*x[0]*x[5]*cos(x[1]));
c3=acos((x[2]*x[2]+x[3]*x[3]-2.385*2.385)/(2*x[2]*x[3])); d1=c1-x[1]-pi/6;
d2=d1-(c2+c3);
aa=x[5]*cos(d1)+x[4]*cos(d2);
bb=x[5]*sin(d1)+x[4]*sin(d2)+9;
f=f+(aa-Dealpos1(n))*(bb-Dealpos2(n));
}
return(f);
}
void strain(double x[],double c[])
{float Lmin=1.772,dlt_L=0.25;
float L;
int n;
for (n=0;n<5;n++)
{L=Lmin+dlt_L*n;
c[1]=pi-acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]));
c[2]=acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]));
c[3]=-acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]))+acos((x[0]*x[0]+x[5]*x[5]-x[1]*x[1])/(2*x[0]*x[5]))+2*pi/3;
c[4]=acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]))+acos((x[0]*x[0]+x[5]*x[5]-x[1]*x[1])/(2*x[0]*x[5]))+pi/3;
c[5]=pi-acos((x[2]*x[2]+x[3]*x[3]-2.385*2.385)/(2*x[2]*x[3]));
c[6]=acos((x[2]*x[2]+x[3]*x[3]-2.385*2.385)/(2*x[2]*x[3]));
c[7]=pi/2-acos((x[1]*x[1]+x[5]*x[5]-x[0]*x[0])/(2*x[1]*x[5]));
c[8]=acos((x[1]*x[1]+x[5]*x[5]-x[0]*x[0])/(2*x[2]*x[5]));
c[9]=x[5]-1-x[4];
}
}
double ldf(double *x)
{int i;
double f,sg;
double c[kkg];
sg=0;
strain(x,c);
for(i=0;i<kkg;i++)
{if(c[i]>0)
sg=sg+r0/c[i];
else
sg=sg-c[i]*1e5;
}
f=objf(x)+sg;
return(f);
}
void ii(double *p,double a[],double b[],double s[],double h0,int n)
{int i;
double *x[3],h,f1,f2,f3;
for(i=0;i<3;i++)
x[i]=(double *)malloc(n*sizeof(double));
h=h0;
for(i=0;i<n;i++)
*(x[0]+i)=*(p+i);
f1=ldf(x[0]);
for(i=0;i<n;i++)
*(x[1]+i)=*(x[0]+i)+h*s[i];
f2=ldf(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)+h*s[i];
f3=ldf(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<3;i++)
free(x[i]);
}
double gold(double a[],double b[],double xx[],int n,double eps1)
{int i;
double f1,f2,*x[2],ff,q,w;
for(i=0;i<2;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=ldf(x[0]);
f2=ldf(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=ldf(x[1]);
}
{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=ldf(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>eps1);
for(i=0;i<n;i++)
xx[i]=0.5*(a[i]+b[i]);
ff=ldf(xx);
for(i=0;i<2;i++)
free(x[i]);
return(ff);
}
double powell(double *p,double x[],double epss,int n)
{int i,j,m;
double *xx[6],*ss,*s,f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
s=(double*)malloc(n*sizeof(double));
ss=(double*)malloc(n*(n+1)*sizeof(double));
for(i=0;i<n;i++)
{for(j=0;j<(n+1);j++)
*(ss+i*(n+1)+j)=0;
*(ss+i*(n+1)+i)=1;
}
for(i=0;i<6;i++)
xx[i]=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
*(xx[0]+i)=*(p+i);
for(;;)
{for(i=0;i<n;i++)
{*(xx[1]+i)=*(xx[0]+i);
x[i]=*(xx[1]+i);
}
f0=f1=ldf(x);
dlt=-1;
for(j=0;j<n;j++)
{for(i=0;i<n;i++)
{*(xx[1]+i)=x[i];
*(s+i)=*(ss+i*(n+1)+j);
}
ii(xx[0],xx[4],xx[5],s,0.1,n);
f=gold(xx[4],xx[5],x,n,1e-5);
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<epss)
{for(i=0;i<6;i++)
free(xx[i]);
free(s);
free(ss);
return(f);
}
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=ldf(x);
f3=fx;
q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
d=0.5*dlt*(f1-f3)*(f1-f3);
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*(n+1)+n+1)=x[i]-*(xx[1]+i);
*(s+i)=*(ss+i*(n+1)+n+1);
}
ii(xx[0],xx[4],xx[5],s,0.1,n);
f=gold(xx[4],xx[5],x,n,1e-5);
for(i=0;i<n;i++)
*(xx[0]+i)=x[i];
for(j=m+1;j<(n+1);j++)
for(i=0;i<n;i++)
*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
}
}
double empf(double *p,double x[],double c,double eps,int n)
{int i,Tm;
double fom,fxo;
fom=1e9;
Tm=0;
do{fxo=powell(p,x,0.01,n);
if(fabs(fom-fxo)>eps)
{fom=fxo;
r0=c*r0;
Tm=Tm+1;
printf("Now it is the %d time of iteration,current values of variables are:\n",Tm);
for(i=0;i<n;i++)
{*(p+i)=x[i];
printf("x[%d}=%f\n",i,x[i]);
}
printf("The value of objective function is:%f\n",fxo);
printf("****************\n");
}
else
return(fxo);
}while(1);
}
void main()
{
int i;
double z[]={2.4,7.5,3.3,1.0,6.5,9.5};
double c,x[6];
c=empf(z,x,0.2,1e-5,6);
printf("\n Output the optimum results:\n");
for(i=0;i<6;i++)
printf("%f\n",x[i]);
printf("The minimum objiective value is %f\n",c);
}
[[it] 本帖最后由 木美子 于 2008-5-26 08:17 编辑 [/it]]
程序代码：