powell 优化算法的过程

发布网友发布时间：2022-11-05 12:56

共1个回答

热心网友时间：2023-10-22 04:46

Powell优化算法是利用仪器测井理建立误差函数（非相关函数），借助Powell方向加速法求出非相关函数达到最小时的解，对于气，水两相流动，从预设的气，水流量初始值出发，沿不同的广向进行搜索，可求出气，水两相流动中可能最大产量。与目前常用的生产测井解释方法相比，文中提出的方法具有精度高，实用性强等优点，在测井曲线有缺陷时，仍有可能得到较好的结果。
powell.c代码如下：
CODE:

＃i nclude "hjfgf.c"
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);
}
}
}
或者%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)<r)
kk=0;
break;
else
ww=0;
m=1;
for i=1:n
gg=ff(y(i,1:n),q)-ff(y(i+1,1:n),q);
if (gg>=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<cha1)
s(n+1)=HJ(y(n+1,1:n),h,a,d(n+1,1:n),q)
xx(k+1,1:2)=y(n+1,1:n)+s(n+1).*d(n+1,1:n)
for j=m+1:n
d(j,1:n)=d(j+1,1:n);
end
k=k+1;
else
xx(k+1,1:n)=y(n+1,1:n);
k=k+1;
end
end
end
powel=y(n+1,1:n)

function w=HJ(x0,h,d,dd,q) %0.618算法
[a,b]=JTF(x0,h,d,dd,q); %调用进退法算法，确定范围
r=0.618;
r1=a+(1-r)*(b-a);
r2=a+r*(b-a);
y1=ff(x0+r1.*dd,q);
y2=ff(x0+r2.*dd,q);
k=1;
while (abs(r1-r2)>=0.1)
if y1<y2
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<y0
h=d*h;
r=r0;
r0=r1;
y0=y1;
else
if k==0;
h=-h;
r=r0;
else
l=0;
break;
end
end
k=k+1;
end
a=min(r,r1);
b=max(r,r1);