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

#include stdio.h#include stdlib.h#include math.hconst int kkg=3;double r0;double f(double x)double ff; ff=pow(x0-8),2)+pow(x1-8),2); return(ff);/*约束条件子程序*/void strain(double x,double g)g0=x0-1; g1=x1-1; g2=11-x0-x1;/*惩罚函数子程序*/double objf(double p)int i; double ff,sg,*g;g=(double *)malloc(kkg*sizeof(double);sg=0;strain(p,g);for(i=0;i0) sg=sg+r0/(*(g+i); elsesg=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 *xx3,h,f1,f2,f3;for (i=0;i3;i+)xxi=(double *)malloc(n*sizeof(double);h=h0;for(i=0;in;i+)*(xx0+i)=x0i;f1=objf(xx0);for(i=0;i=f1)h=-h0;for(i=0;in;i+)*(xx2+i)=*(xx0+i);f3=f1; for(i=0;in;i+)*(xx0+i)=*(xx1+i); *(xx1+i)=*(xx2+i);f1=f2;f2=f3;for(;)h=2.0*h;for(i=0;in;i+)*(xx2+i)=*(xx1+i)+h*si;f3=objf(xx2);if(f2f3)break;elsefor(i=0;in;i+)*(xx0+i)=*(xx1+i); *(xx1+i)=*(xx2+i);f1=f2;f2=f3;if(h0)for(i=0;in;i+)ai=*(xx2+i); bi=*(xx0+i);elsefor(i=0;in;i+)ai=*(xx0+i); bi=*(xx2+i);for(i=0;i3;i+)free(xxi);/*黄金分割程序*/double gold (double a,double b,double eps,int n,double x)int i;double f1,f2,*xx2,ff,q,w;for(i=0;i2;i+)xxi=(double *)malloc(n*sizeof(double);for(i=0;if2)for(i=0;in;i+)bi=*(xx0+i); *(xx0+i)=*(xx1+i);f1=f2;for(i=0;in;i+)*(xx1+i)=ai+0.382*(bi-ai);f2=objf(xx1);elsefor(i=0;in;i+)ai=*(xx1+i); *(xx1+i)=*(xx0+i);f2=f1;for(i=0;in;i+)*(xx0+i)=ai+0.618*(bi-ai);f1=objf(xx0);q=0;for(i=0;ieps);for(i=0;in;i+)xi=0.5*(ai+bi);ff=objf(x);for(i=0;i2;i+)free(xxi);return(ff);double oneoptim(double x0,double s,double h0,double epsg,int n,double x)double ff,*a,*b; 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);/*powell子程序*/double powell(double p,double h0,double eps,double epsg,int n,double x)int i,j,m; double *xx4,*s,*ss; 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;i4;i+) xxi=(double*)malloc(n*sizeof(double); for(i=0;in;i+) *(xx0+i)=pi; for(;) for(i=0;in;i+) *(xx1+i)=*(xx0+i); xi=*(xx1+i); f0=f1=objf(x); dlt=-1; for(j=0;jn;j+) for(i=0;idlt) dlt=df; m=j; sdx=0; for(i=0;in;i+) sdx=sdx+fabs(xi-(*(xx1+i); if(sdxeps) for(i=0;i4;i+) free(xxi); return(f); for(i=0;in;i+) *(xx2+i)=xi; f2=f; for(i=0;in;i+) *(xx3+i)=2*(*(xx2+i)-(*(xx1+i); xi=*(xx3+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(f3f1)|(qd) if(f2=f3) for(i=0;in;i+) *(xx0+i)=*(xx2+i); else for(i=0;in;i+) *(xx0+i)=*(xx3+i); else for(i=0;in;i+) *(ss+(i+1)*(n+1)=xi-(*(xx1+i); *(s+i)=*(ss+(i+1)*(n+1); f=oneoptim(xx0,s,h0,epsg,n,x); for(i=0;in;i+) *(xx0+i)=xi; for(j=m+1;j=