论坛 › 算法及编程语言 › 内点法优化程序

donglei06 发表于 2007-5-18 14:50

内点法优化程序

内点法优化程序
#include<math.h>
#include<stdlib.h>
#include<stdio.h>
const int kkg=8;
double r0;
void main()
{double f(double x[]);
double objf(double p[]);
void jtf(double x0[],double h0, double s[],int n, double a[],double b[]);
double gold(double a[],double b[],double eps,int n, double xx[]);
double oneoptim(double x0[],double s[],double h0, double epsg,int n, double x[]);
double powell(double p[],double h0, double eps, double epsg,int n, double x[]);
int i;
double p[]={21,35,108,16};
double fom,fxo,c,x;
c=0.5;
r0=45800;
fom=10000;
do
{
fxo=powell(p,0.1,1500,650,4,x);
if(fabs(fom-fxo)>1800)
{fom=fxo;
r0=c*r0;
for(i=0;i<4;i++)
*(p+i)=x;
}
else
{printf("x=%f,x=%f,x=%f,x=%f\n",x,x,x,x);
return;
}
}while(1);
}
double f(double x[])
{ double ff;
ff=0.8322*x*(3.7889*x*x+5.8174*x+3.9422)*pow((0.9715*x*x-0.5*x),2)
/(cos(0.0087*x)*x*x);
return(ff);
}
void strain(double x[], double g[])
{g=4.842e-6*x*x*x*x-1.0;
g=4.86e-4*x*x*sqrt(x)-1.0;
g=x-12;
g=x-30;
g=45-x;
g=x-1.5;
g=x-0.2429*x*x;
g=0.3238*x*x-x;
}
double objf(double p[])
{int i;
double ff,sg,*g;
g=(double*)malloc(kkg*sizeof(double));
sg=0;
for(i=0;i<kkg;i++)
{if(*(g+i)>0)
sg=sg+r0/(*(g+i));
else
sg=sg+r0*(1e+10);
}
free(g);
ff=f(p)+sg;
return(ff);
}
void jtf(double x0[],double h0, double s[],int n, double a[],double b[])
{int i;
double *x,h,f1,f2,f3;
for(i=0;i<3;i++)
x=(double*)malloc(n*sizeof(double));
h=h0;
for(i=0;i<n;i++)
   *(x+i)=x0;
   f1=objf(x);
   for(i=0;i<n;i++)
   *(x+i)=*(x+i)+h*s;
   f2=objf(x);
if(f2>=f1)
   {h=-h0;
   for(i=0;i<n;i++)
   *(x+i)=*(x+i);
   f3=f1;
   for(i=0;i<n;i++)
      {*(x+i)=*(x+i);
      *(x+i)=*(x+i);
      }
   f1=f2;
   f2=f3;
   }
for(;;)
{h=2*h;
   for(i=0;i<n;i++)
   *(x+i)=*(x+i)+h*s;
   f3=objf(x);
   if(f2<f3)
   break;
else
{for(i=0;i<n;i++)
{*(x+i)=*(x+i);
   *(x+i)=*(x+i);
}
f1=f2;
f2=f3;
}
}
if(h<0)
for(i=0;i<n;i++)
{a=*(x+i);
   b=*(x+i);
}
else
for(i=0;i<n;i++)
{a=*(x+i);
b=*(x+i);
}
for(i=0;i<3;i++)
free(x);
}
double gold(double a[],double b[],double eps,int n, double xx[])
{int i;
double f1,f2,*x,ff,q,w;
for(i=0;i<2;i++)
x=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{*x=a+0.618*(b-a);
   *(x+i)=a+0.382*(b-a);
}
f1=objf(x);
f2=objf(x);
do
{if(f1>f2)
   { for(i=0;i<n;i++)
   {b=*(x+i);
      *(x+i)=*(x+i);
}
f1=f2;
for(i=0;i<n;i++)
   *(x+i)=a+0.382*(b-a);
   f2=objf(x);
}
else
{ for(i=0;i<n;i++)
   {a=*(x+i);
   *(x+i)= *(x+i);
    }
f2=f1;
for(i=0;i<n;i++)
*x=a+0.618*(b-a);
f1=objf(x);
}
q=0;
for(i=0;i<n;i++)
q=q+(b-a)*(b-a);
w=sqrt(q);
}while(w>eps);
for(i=0;i<n;i++)
xx=0.5*(a+b);
ff=objf(xx);
for(i=0;i<2;i++)
free(x);
return(ff);
}
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
{double *a,*b,ff;
a=(double*)malloc(n*sizeof(double));
b=(double*)malloc(n*sizeof(double));
jtf(x0,h0,s,n,a,b);
ff=gold(a,b,epsg,n,x);
free(a);
free(b);
return(ff);
}
double powell(double p[],double h0,double eps,double epsg,int n,double x[])
{int i,j,m;
double *xx,*ss,*s;
double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
ss=(double*)malloc(n*(n+1)*sizeof(double));
s=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{for(j=0;j<=n;j++)
*(ss+i*(n+1)+j)=0;
*(ss+i*(n+1)+i)=1;
}
for(i=0;i<4;i++)
xx=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
*(xx+i)=p;
for(;;)
{for(i=0;i<n;i++)
{*(xx+i)=*(xx+i);
x=*(xx+i);
}
f0=f1=objf(x);
dlt=-1;
for(j=0;j<n;j++)
{for(i=0;i<n;i++)
   {*(xx+i)=x;
    *(s+i)=*(ss+i*(n+1)+j);
   }
f=oneoptim(xx,s,h0,epsg,n,x);
df=f0-f;
if(df>dlt)
{dlt=df;
m=j;
}
}
sdx=0;
for(i=0;i<n;i++)
sdx=sdx+fabs(x-(*(xx+i)));
if(sdx<eps)
{free(ss);
   free(s);
   for(i=0;i<4;i++)
    free(xx);
   return(f);
}
for(i=0;i<n;i++)
*(xx+i)=x;
f2=f;
for(i=0;i<n;i++)
{*(xx+i)=2*(*(xx+i)-(*(xx+i)));
   x=*(xx+i);
}
fx=objf(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+i)=*(xx+i);
      else
      for(i=0;i<n;i++)
         *(xx+i)=*(xx+i);
   }
else
{for(i=0;i<n;i++)
   {*(ss+(i+1)*(n+1))=x-(*(xx+i));
    *(s+i)=*(ss+(i+1)*(n+1));
   }
f=oneoptim(xx,s,h0,epsg,n,x);
for(i=0;i<n;i++)
   *(xx+i)=x;
for(j=m+1;j<=n;j++)
for(i=0;i<n;i++)
   *(ss+i*(n+1)+j-i)=*(ss+i*(n+1)+j);
}
}
}
运行结果是:Floating point error:Over flow

[ 本帖最后由 风花雪月 于 2007-6-15 11:15 编辑 ]

hnsqzrj 发表于 2007-6-8 20:37

能运行么

总有一天 发表于 2007-6-15 01:10

对不起，我觉得还挺长的，建议楼主以后加上注释

风花雪月 发表于 2007-6-15 11:19

自己好好找找是什么地方溢出，我运行了一下，好像需要运行很长时间

估计这样的问题，别人是没有那么多时间帮你调试的

查看完整版本: 内点法优化程序