复合单纯型法C语言程序

#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) //定义gau函数 { int i; constraint(x,g); for(i=0;i<kg;i++) { if(g[i]<0) goto s333; /*gau=1则为可行点，否则为非可行点*/ } return 1; s333:return 0; } void xcent(int n,int k,int ll,int lh,float x0[],float xcom[][4]) { int i,l; float xs; for(i=0;i<n;i++) { xs=0; for(i=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[][4],float fxk[]) { int l,lp,lp1,i; float w; for(l=0;l<k-1;l++) for(lp=0;lp<k-1;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]=x[i]; } } } } void complex(int n,int k,int kg,float ep,float x[],float bl[],float bu[],float xcom[][4],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); } float objfx(float x[]) { return x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60; //目标函数 } void constraint(float x[],float g[]) //约束条件 { g[0]=x[0]; g[1]=x[1]; g[2]=6.0-x[0]; g[3]=8.0-x[1]; g[4]=11.0-x[0]-x[1]; } void main() { float x[2],bl[2],bu[2],xcom[2][4],f; int n=2,nf=0,ng=0; bl[0]=0.0; bl[1]=0.0; bu[0]=9.0; bu[1]=9.0; system("cls"); complex(n,4,5,1.0e-5,x,bl,bu,xcom,&f,&nf,&ng); printf("\n x[0]=%f,x[1]=%f",x[0],x[1]); //输出最优解 printf("\n minF(x)=%f",f); //输出最优值 printf("\n nf=%d,ng=%d",nf,ng); getchar(); }