内点惩罚函数法子程序.txt

内点惩罚函数法,机械优化设计中的优化方法.

void gradient(double x[],double g[],int n)
{int i;
  double af,f1,f2,dltx=0.000001;
  for(i=0;i<n;i++)
  g[i]=0;
  f1=objf(x);
  for(i=0;i<n;i++)
  {af=*(x+i);
   *(x+i)=af+dltx;
   f2=objf(x);
   g[i]=(f2-f1)/dltx;
   *(x+i)=af;
  }
}
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 dfpopt(double xx[],double h0,double eps,double epsg,int n)
{int i,j,k;
  double ae,zcc;
  double *s,*x,*ay[2],*df[2],*zd[2],*zc[2],*zh[2];
  s=(double *)malloc(n*sizeof (double ));
  x=(double *)malloc(n*sizeof (double ));
  for(i=0;i<2;i++)
  {ay[i]=(double*)malloc(n*sizeof(double));
   df[i]=(double*)malloc(n*sizeof(double));
   zd[i]=(double*)malloc(n*sizeof(double));
   zc[i]=(double*)malloc(n*sizeof(double));
   zh[i]=(double*)malloc(n*n*sizeof(double));
  }
  for(i=0;i<n;i++)
  *(ay[0]+i)=xx[i];
  L1:
   k=0;
   for (i=0;i<n;i++)
   for (j=0;j<n;j++)
    {*(zh[0]+i*n+j)=0;
     if(i==j)
      *(zh[0]+i*n+j)=1;
    }
  for(i=0;i<n;i++)
  {*(x+i)=*(ay[0]+i);
   *(ay[1]+i)=*(ay[0]+i);
  }
  gradient(x,df[0],n);
  for(i=0;i<n;i++)
  {*(s+i)=0;
   *(df[1]+i)=*(df[0]+i);
   for(j=0;j<n;j++)
   *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
  }
  L2:
  oneoptim(x,s,h0,epsg,n,ay[0]);
  for(i=0;i<n;i++)
   *(x+i)=*(ay[0]+i);
  gradient(x,df[0],n);
  ae=0;
  for(i=0;i<n;i++)
   ae=ae+(*(df[0]+i))*(*(df[0]+i));
  if(ae<=eps)
   {for(i=0;i<n;i++)
    *(xx+i)=*(x+i);
    free(s);
    free(x);
    for(i=0;i<2;i++)
    {free(ay[i]);
     free(df[i]);
     free(zd[i]);
     free(zc[i]);
     free(zh[i]);
    }
    return(objf(xx));
   }
   if(k==n) goto L1;
   zcc=0;
   for(i=0;i<n;i++)
   {*(zd[0]+i)=*(ay[0]+i)-(*(ay[1]+i));
    *(zd[1]+i)=*(df[0]+i)-(*(df[1]+i));
    *(df[1]+i)=*(df[0]+i);
    zcc=zcc+(*(zd[0]+i))*(*(zd[1]+i));
   }
   for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    *(zh[1]+i*n+j)=*(zd[0]+i)*(*(zd[0]+j))/zcc;
   for(i=0;i<n;i++)
    { *(zc[0]+i)=0;
      for(j=0;j<n;j++)
       *(zc[0]+i)=*(zc[0]+i)+(*(zd[1]+j))*(*(zh[0]+j*n+i));
    }
     zcc=0;
    for(i=0;i<n;i++)
     zcc=zcc+(*(zc[0]+i))*(*(zd[1]+i));
     for(i=0;i<n;i++)
     {*(zc[0]+i)=0;
      *(zc[1]+i)=0;
      for (j=0;j<n;j++)
       {*(zc[0]+i)=*(zc[0]+i)+(*(zh[0]+i*n+j))*(*(zd[1]+j));
        *(zc[1]+i)=*(zc[1]+i)+(*(zh[1]+j))*(*(zh[0]+j*n+i));
       }
     }
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    *(zh[0]+i*n+j)=*(zh[0]+i*n+j)+(*(zh[1]+i*n+j))-(*(zc[0]+i))*(*(zc[1]+j))/zcc;
    for(i=0;i<n;i++)
     {*(s+i)=0;
      for(j=0;j<n;j++)
       *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
     }
     k=k+1;
     goto L2;
}



#include "stdlib.h"
#include "math.h"
const int kkg=3;
double r0;
double f(double x[])
{double ff;
ff=pow((x[0]-8),2)+pow((x[1])-8,2);
return(ff);
}
void strain(double  x[],double g[])
{   g[0]=x[0]-1;
    g[1]=x[1]-1;
    g[2]=11-x[0]-x[1];
}
double objf(double p[])
{int i;
double ff,sg,*g;
g=(double *)malloc(kkg *sizeof(double));
sg=0;
strain(p,g);
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);
}



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[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=objf(x[0]);
  f2=objf(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=objf(x[1]);
     }
    else
     { 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=objf(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>eps);
  for(i=0;i<n;i++)
   xx[i]=0.5*(a[i]+b[i]);
  ff=objf(xx);
  for(i=0;i<2;i++)
  free(x[i]);
  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[i]=(double *)malloc(n*sizeof(double));
  h=h0;
for(i=0;i<n;i++)
  *(x[0]+i)=x0[i];
f1=objf(x[0]);
for(i=0;i<n;i++)
  *(x[1]+i)=*(x[0]+i)+h*s[i];
  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[i];
   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[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 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[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=objf(x);
   dlt=-1;
   for(j=0;j<n;j++)
    {for(i=0;i<n;i++)
      {*(xx[0]+i)=x[i];
       *(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[i]-(*(xx[1]+i)));
    if(sdx<eps)
     {free(ss);
      free(s);
      for(i=0;i<4;i++)
      free(xx[i]);
      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=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[i]-(*(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[i];
   for(j=m+1;j<=n;j++)
   for(i=0;i<n;i++)
    *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
    }
  }
}



main()
{int i;
double p[]={3,4};
double fom,fxo,x[2];
double c=0.1;
double  r0=120;
fom=100;
do
  { fxo=dfpopt(p,0.2,0.0001,0.000001,2);
    if(fabs(fom-fxo)>0.001)
    {fom=fxo;
     r0=c*r0;
     for(i=0;i<2;i++)
     *(p+i)=x[i];
    }
   else
    {printf("\nx[0]=%f,x[1]=%f,ff=%f",x[0],x[1],fxo);
     return; }
    }while(1);
    getch();
}