1.c

用VC编写的内点罚函数法

/* 以下为复合形法法子程序 */

#include "stdlib.h"

#include "math.h"

#include "stdio.h"

float objfx(float x[]);

void constraint(float x[],float g[]);
int gau(float x[],float g[],int kg)

{

       int i;

       constraint(x,g);

       for(i=0;i<kg;i++)

       {

         if(g[i]<0)

         goto s333;

       }

       return 1;

s333:return 0;

}
void xcent(int n,int k,int ll,int lh,float x0[],float xcom[][100])

{

       int i,l;

       float xs;

       for(i=0;i<n;i++)

       {

         xs=0;

         for(l=0;l<ll;l++)

         {

    if(l!=lh)

    xs=xs+xcom[i][l];

         }

         if(lh>-1)

    x0[i]=xs/(ll-1);

         else

    x0[i]=xs/ll;

       }

}
void fxse(int n,int k,float x[],float xcom[][100],float fxk[])

{

       int l,lp,lp1,i;

       float w;

       for(l=0;l<k-1;l++)

       for(lp=0;lp<k-l;lp++)

       {

         lp1=lp+1;

         if(fxk[lp]<=fxk[lp1])

         {

    w=fxk[lp];

    fxk[lp]=fxk[lp1];

    fxk[lp1]=w;

    for(i=0;i<n;i++)

    {

      x[i]=xcom[i][lp];

      xcom[i][lp]=xcom[i][lp1];

      xcom[i][lp1]=x[i];

    }

         }

       }

}


void complex(int n,int k,int kg,float ep,float x[],float bl[],float bu[],

        float xcom[][100],float *f,int *nf,int *ng)

{

       int i,iw,l,ll,lh,it;

       float fx,fx0,sdx,fxh,fxr,alp;

       float *x0=(float*)calloc(n,sizeof(float));

       float *xh=(float*)calloc(n,sizeof(float));

       float *xr=(float*)calloc(n,sizeof(float));

       float *fxk=(float*)calloc(k,sizeof(float));

       float *g=(float*)calloc(kg,sizeof(float));

s5:    for(i=0;i<n;i++)

       x[i]=bl[i]+rand()/40000.0*(bu[i]-bl[i]);

       iw=gau(x,g,kg);

       *ng=*ng+1;

       if(iw==0)goto s5;

       for(i=0;i<n;i++)

       xcom[i][0]=x[i];

       for(l=1;l<k;l++)

       for(i=0;i<n;i++)

       xcom[i][l]=bl[i]+rand()/50000.0*(bu[i]-bl[i]);

       lh=-1;

       for(ll=1;ll<k;ll++)

       {

         xcent(n,k,ll,lh,x0,xcom);

         iw=gau(x0,g,kg);

         *ng=*ng+1;

         if(iw==0)goto s5;

         for(i=0;i<n;i++)

         x[i]=xcom[i][ll+1];

s24:     iw=gau(x,g,kg);

         *ng=*ng+1;

         if(iw==0)

         {

    for(i=0;i<n;i++)

    x[i]=x0[i]+0.5*(x[i]-x0[i]);

    goto s24;

         }

         else

         {

    for(i=0;i<n;i++)

    xcom[i][ll+1]=x[i];

         }

       }

       for(l=0;l<k;l++)

       {

         for(i=0;i<n;i++)

         x[i]=xcom[i][l];

         fx=objfx(x);

         *nf=*nf+1;

         fxk[l]=fx;

       }

       it=0;

s14: it=it+1;

       printf("\n\n    +++ ITERATION COMPUTE +++");

       printf("\n       ITER    = %d",it);

       lh=-1;

       xcent(n,k,k,lh,x0,xcom);

       fx0=objfx(x0);

       *nf=*nf+1;

       iw=gau(x0,g,kg);

       *ng=*ng+1;

       printf("\n       Fmid    = %f",fx0);

       for(i=0;i<n;i++)

         printf("\n       X(%d)mid=%f",i,x0[i]);

       for(i=0;i<kg;i++)

         printf("\n       G(%d)mid=%f",i,g[i]);

       sdx=0;

       for(l=0;l<k;l++)

         sdx=sdx+(fx0-fxk[l])*(fx0-fxk[l]);

       sdx=sqrt(sdx/(float)k);

       if(sdx<ep)goto s38;

       fxse(n,k,x,xcom,fxk);

       lh=0;

s22: fxh=fxk[lh];

       for(i=0;i<n;i++)

         xh[i]=xcom[i][lh];

       xcent(n,k,k,lh,x0,xcom);

       iw=gau(x0,g,kg);

       *ng=*ng+1;

       if(iw==0)goto s36;

       alp=1.3;

s12: for(i=0;i<n;i++)

       xr[i]=x0[i]+alp*(x0[i]-xh[i]);

       iw=gau(xr,g,kg);

       *ng=*ng+1;

       if(iw==0)

       {

         alp=alp*0.5;

         goto s12;

       }

       fxr=objfx(xr);

       *nf=*nf+1;

       if(fxr>=fxh)

       {

         if(alp>1.0e-4)

         {

    alp=alp*0.5;

    goto s12;

         }

         lh=lh+1;

         if(lh<3)goto s22;

       }

       for(i=0;i<n;i++)

       xcom[i][lh]=xr[i];

       fxk[lh]=fxr;

       goto s14;

s36: for(i=0;i<n;i++)

       {

         bl[i]=xcom[i][k];

         bu[i]=x0[i];

       }

       goto s5;

s38: for(i=0;i<n;i++)

       x[i]=x0[i];

       *f=objfx(x);

       *nf=*nf+1;

       free(x0);

       free(xh);

       free(xr);

       free(g);

       free(fxk);

}