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=0i<ni++)
{for(j=0j<=nj++)
*(ss+i*(n+1)+j)=0
*(ss+i*(n+1)+i)=1
}
for(i=0i<4i++)
xx[i]=(double *)malloc(n*sizeof(double))
for(i=0i<ni++)
*(xx[0]+i)=p[i]
for()
{for(i=0i<ni++)
{*(xx[1]+i)=*(xx[0]+i)
x[i]=*(xx[1]+i)
}
f0=f1=objf(x)
dlt=-1
for(j=0j<nj++)
{for(i=0i<ni++)
{*(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=0i<ni++)
sdx=sdx+fabs(x[i]-(*(xx[1]+i)))
if(sdx<eps)
{free(ss)
free(s)
for(i=0i<4i++)
free(xx[i])
return(f)
}
for(i=0i<ni++)
*(xx[2]+i)=x[i]
f2=f
for(i=0i<ni++)
{*(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=0i<ni++)
*(xx[0]+i)=*(xx[2]+i)
else
for(i=0i<ni++)
*(xx[0]+i)=*(xx[3]+i)
}
else
{for(i=0i<ni++)
{*(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=0i<ni++)
*(xx[0]+i)=x[i]
for(j=m+1j<=nj++)
for(i=0i<ni++)
*(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)
鲍威尔基本算法的问题在于，可能发生退化问题，具体而言就是可能在某一环迭代中出现基本方向组线性相关的情况，这种情况下按新方向替代第一个方向的方法进行替换，就会导致搜索在降维的空间中进行，无法得到原本n维空间的函数极小值，计算将失败。
而改进的方法和原来方法本质区别在于替换方向的规则不同。改进的方法，能够保证每轮迭代中搜索方向都线性无关，而且随着迭代的延续，共轭的程度会逐渐增加。
具体展开比较复杂，简单来说就是每次产生了新生方向，都要判断一下这个方向好不好，如果不好就不换进来；如果觉得这个方向好，就看一下旧方向中哪个函数下降量最大，把这个下降量最大的方向替换掉。
