[计算机]Powell优化算法是利用仪器测井理建立误差函数.doc

Powell优化算法是利用仪器测井理建立误差函数（非相关函数），借助Powell方向加速法求出非相关函数达到最小时的解，对于气，水两相流动，从预设的气，水流量初始值出发，沿不同的广向进行搜索，可求出气，水两相流动中可能最大产量。与目前常用的生产测井解释方法相比，文中提出的方法具有精度高，实用性强等优点，在测井曲线有缺陷时，仍有可能得到较好的结果。powell.c代码如下： CODE:i nclude hjfgf.cdouble 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 *xx4,*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;in;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) free(ss); free(s); 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=n;j+) for(i=0;in;i+) *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j); 或者%powell算法，用于寻找无约束最优值点 function powel=powell(x0,n,q,r,h,a) d=eye(n); %n个线性无关的初始搜索方向 k=1; kk=1; xx(1,1:n)=x0; while (kk) y(1,1:n)=xx(k,1:n); for j=1:n s(j)=HJ(y(j,1:n),d(j,1:n),q); %调用0.618算法 y(j+1,1:n)=y(j,1:n)+s(j).*d(j,1:n); end d(n+1,1:n)=y(n+1,1:n)-y(1,1:n); if (norm(d(n+1,1:n),2)=ww) m=i; end end cha=ff(y(1,1:n),q)-2*ff(y(n+1,1:n),q)+ff(2*y(n+1,1:n)-y(1,1:n),q); cha1=2*(ff(y(m,1:n),q)-ff(y(m+1,1:n),q); if (cha=0.1) if y1y2 b=r2; r2=r1; y2=y1; r1=a+(1-r)*(b-a); y1=ff(x0+r1.*dd,q); else a=r1; r1=r2; y1=y2; r2=a+r*(b-a); y2=ff(x0+r2.*dd,q); end end w=(r1+r2)/2%进退法function a,b=JTF(x0,h,d,dd,q) r0=0; y0=ff(x0+r0.*dd,q); k=0; l=1; while (l) r1=r0+h; y1=ff(x0+r1.*dd,q); if y1