新建 文本文档 (2).txt

程序在vc++6.0下面通过

内点法优化程序
#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[4];
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[0]=%f,x[1]=%f,x[2]=%f,x[3]=%f\n",x[0],x[1],x[2],x[3]);
return;
}
}while(1);
}
double f(double x[])
{ double ff;
ff=0.8322*x[2]*(3.7889*x[0]*x[0]+5.8174*x[0]+3.9422)*pow((0.9715*x[3]*x[0]-0.5*x[2]),2)
/(cos(0.0087*x[1])*x[0]*x[0]);
return(ff);
}
void strain(double x[], double g[])
{g[0]=4.842e-6*x[0]*x[2]*x[3]*x[3]-1.0;
g[1]=4.86e-4*x[0]*x[3]*sqrt(x[2])-1.0;
g[2]=x[0]-12;
g[3]=x[1]-30;
g[4]=45-x[1];
g[5]=x[3]-1.5;
g[6]=x[2]-0.2429*x[3]*x[0];
g[7]=0.3238*x[3]*x[0]-x[2];
}
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[3],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[0]+i)=x0;
   f1=objf(x[0]);
   for(i=0;i<n;i++)
   *(x[1]+i)=*(x[0]+i)+h*s;
   f2=objf(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;
   f3=objf(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=*(x[2]+i);
   b=*(x[0]+i);
  }
else
for(i=0;i<n;i++)
{a=*(x[0]+i);
  b=*(x[2]+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[2],ff,q,w;
  for(i=0;i<2;i++)
  x=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
  {*x[0]=a+0.618*(b-a);
   *(x[1]+i)=a+0.382*(b-a);
  }
f1=objf(x[0]);
f2=objf(x[1]);
do
{if(f1>f2)
   { for(i=0;i<n;i++)
     {b=*(x[0]+i);
      *(x[0]+i)=*(x[1]+i);
  }
f1=f2;
for(i=0;i<n;i++)
   *(x[1]+i)=a+0.382*(b-a);
   f2=objf(x[1]);
}
else
{ for(i=0;i<n;i++)
   {a=*(x[1]+i);
     *(x[1]+i)= *(x[0]+i);
    }
  f2=f1;
for(i=0;i<n;i++)
*x[0]=a+0.618*(b-a);
f1=objf(x[0]);
}
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[4],*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[0]+i)=p;
for(;;)
{for(i=0;i<n;i++)
{*(xx[1]+i)=*(xx[0]+i);
  x=*(xx[1]+i);
  }
f0=f1=objf(x);
dlt=-1;
for(j=0;j<n;j++)
  {for(i=0;i<n;i++)
   {*(xx[0]+i)=x;
    *(s+i)=*(ss+i*(n+1)+j);
   }
f=oneoptim(xx[0],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[1]+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[2]+i)=x;
f2=f;
for(i=0;i<n;i++)
  {*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
   x=*(xx[3]+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[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+1)*(n+1))=x-(*(xx[1]+i));
    *(s+i)=*(ss+(i+1)*(n+1));
   }
  f=oneoptim(xx[0],s,h0,epsg,n,x);
  for(i=0;i<n;i++)
   *(xx[0]+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);
}
}
}