powell法matlab代码

1/目标函数程序清单
function m=f(x1,x2)
m=x1^2+2*x2^2-4*x1-2*x1*x2

2/关于α的目标函数
function m=y(x1,x2,d1,d2,alpha)
m=(x1+alpha*d1)^2+(x2+alpha*d2)^2-(x1+alpha*d1)*(x2+alpha*d2)-10*(x1+alpha*d1)-4*(x2+
alpha*d2)+60;

3\
function [a,b]=section(x1,x2,d1,d2)
x11=x1; %给出起始点坐标x1，x2和搜索方向d1，d2
x22=x2;
d11=d1;
d22=d2;
h0=1; %初始化
h=h0;alpha1=0;
y1=y(x11,x22,d11,d22,alpha1); %代入α1求解y1
alpha2=h;y2=y(x11,x22,d11,d22,alpha2);
t=0;
if y2>y1
h=-h;alpha3=alpha1;y3=y1;t=1; %如果y2>y1，则改变搜索方向
end
while(1)
if t==1
alpha1=alpha2;y1=y2; %实现交换
alpha2=alpha3;y2=y3;
else t=1;
end
alpha3=alpha2+h;y3=y(x11,x22,d11,d22,alpha3);
if y3h=2*h; %改变搜索步长，将其加倍
else
break;
end
end
if alpha1>alpha3
tem=alpha1;alpha1=alpha3;alpha3=tem; %比较大小，保证输出区间为【a，b】
a=alpha1;b=alpha3;
else a=alpha1;b=alpha3;
end

4\
function alpha=ALPHA(x1,x2,d1,d2,A,B)
x11=x1;x22=x2; %给出起始点坐标x1，x2和搜索方向d1，d2
d11=d1;d22=d2;
a=A;b=B; %获取区间
epsilon=0.000001; %初始化，给定进度
r=0.618;
alpha1=b-r*(b-a);
y1=y(x11,x22,d11,d22,alpha1); %代入α1求解y1
alpha2=a+r*(b-a);
y2=y(x11,x22,d11,d22,alpha2);
while(1)
if y1>=y2 %根据区间消去法原理缩短搜索空间
a=alpha1;alpha1=alpha2;
y1=y2;
alpha2=a+r*(b-a);
y2=y(x11,x22,d11,d22,alpha2);
else
b=alpha2;alpha2=alpha1;
y2=y1;
alpha1=b-r*(b-a);
y1=y(x11,x22,d11,d22,alpha1);
end
if abs(b-a)abs(y2-y1)break; %满足进度要求则退出循环迭代过程
end
end
alpha=0.5*(a+b); %返回值

5\主程序
clear all
clc
k=0;n=2;
x=[0;0;];ff(1)=f(x(1),x(2)); %初始化
epsilon=0.00001;
d=[1;
0;
0;
1;];
while(1)
x00=[x(1);x(2);];
for i=1:n
[a(i),b(i)]=section(x(2*i-1),x(2*i),d(2*i-1),d(2*i)); %调用section（）函数求解一元函数有y（α）最小值时的区间
alpha(i)=ALPHA(x(2*i-1),x(2*i),d(2*i-1),d(2*i),a(i),b(i)); %调用ALPHA（）函数求解y（α）最小值点α*
x(2*i+1)=x(2*i-1)+alpha(i)*d(2*i-1); %沿di方向搜索
x(2*i+2)=x(2*i)+alpha(i)*d(2*i);
ff(i+1)=f(x(2*i+1),x(2*i+2)); %搜索到的点对应的函数值
end

for i=1:n
Delta(i)=ff(i)-ff(i+1);
end
delta=max(Delta); %求出函数值之差的最大值
for i=1:n %寻找函数值之差最大值的下标
if delta==Delta(i)
m=i;
break;
end
end
d(2*n+1)=x(2*n+1)-x(1); %求出反射点搜索方向dn+1
d(2*n+2)=x(2*n+2)-x(2);
x(2*n+3)=2*x(2*n+1)-x(1); %搜索到反射点Xn+1
x(2*n+4)=2*x(2*n+2)-x(2);
ff(n+2)=f(x(2*n+3),x(2*n+4)); %反射点所对应的函数

值
f0=ff(1);f2=ff(n+1);f3=ff(n+2);
k=k+1; %记录迭代次数
R(k,:)=[k,x',d',ff]; %保存迭代过程的中间运算结果

if f3(f0-2*f2+f3)*(f0-f2-delta)^2<0.5*delta*(f0-f3)^2 %判断是否需要对原方向组进行替换
[a(n+1),b(n+1)]=section(x(2*n+1),x(2*n+2),d(2*n+1),d(2*n+2));
alpha(n+1)=ALPHA(x(2*n+1),x(2*n+2),d(2*n+1),d(2*n+2),a(n+1),b(n+1));
x(1)=x(2*n+1)+alpha(n+1)*d(2*n+1); %沿反射方向进行搜索，将搜索结果作为下一轮迭代的始点
x(2)=x(2*n+2)+alpha(n+1)*d(2*n+2);
for i=m:n %根据函数值之差最大值的下标值m，对原方向组进行替换
d(2*i-1)=d(2*i+1);
d(2*i)=d(2*i+2);
end
else
if f2x(1)=x(2*n+1);
x(2)=x(2*n+2);
else
x(1)=x(2*n+3);
x(2)=x(2*n+4);
end
end
RR(k,:)=alpha; %保存迭代过程的中间运算结果
ff(1)=f(x(1),x(2)); %计算下一轮迭代过程需要的f0值
if
(((x(2*n+1)-x00(1))^2+(x(2*n+2)-x00(2))^2)^(1/2))break; %满足进度要求则退出循环迭代过程
end
end
xx=[x(1);x(2)]
fmin=f(x(1),x(2)) %显示最小值及其所对应的坐标