用POWELL法求极小值：Rosenbrock函数

函数：

一、题目分析：

       在数学最优化中，Rosenbrock 函数是一个用来测试最优化算法性能的非凸函数，由 Howard Harry Rosenbrock 在 1960 年提出 。也称为 Rosenbrock 山谷或 Rosenbrock 香 蕉函数，也简称为香蕉函数。Rosenbrock 函数的定义如下：

       Rosenbrock函数的每个等高线大致呈抛物线形，其全域最小值也位在抛物线形的山谷中（香蕉型山谷）。

 函数图像

 二、基本方法原理

鲍威尔(Powell)法是直接利用函数值来构造共辄方向的一种共辄方向法。其基本思想是:在不用导数的前提下，在迭代中逐次构造共辄方向。鲍威尔的基本算法是:在每一轮迭代中总有一个始点(第一轮的始点是任选的初始点)和n个线性独立的搜索方向。从始点出发顺次沿n个方向作一维搜索得一终点，由始点和终点决定了一个新的搜索方向。用这个方向替换原来n个方向中的一个，于是形成新的搜索方向组。替换的原则是去掉原方向组的第一个方向而将新方向排在原方向的最后。此外规定，从这一轮的搜索终点出发沿新的搜索方向作一维搜索而得到的极小点，作为下一轮迭代的始点。这样就形成算法的循环。因为这种方法在迭代中逐次生成共辄方向，而共辄方向是较好的搜索方向，所以鲍威尔法又称方向加速法。

针对二维情况,鲍威尔法的寻优过程如图所示。

 图 二维正定函数得Powell法

程序

主程序powell.m

clear all
clc
h0=0.1;
x0=[0,0];
e=[1,0;0,1];
eps = 1e-6;
n=2;
for k=1:100000
 for i=1:n 
 x01=x0(1);x02=x0(2);d=e(i,:); 
 F0=ff(x0(1),x0(2));
 opt_step=goldsection(x01,x02,d,h0);
 x(i,:)=x0+opt_step*d; 
 F(i)=ff(x(i,1),x(i,2));
 D(i)=F(i)-F0;
 end
 for i=2:n 
 x01=x(i-1,1);x02=x(i-1,2);d=e(i,:);
 opt_step=goldsection(x01,x02,d,h0);
 x(i,:)=x(i-1,:)+opt_step*d; 
 F(i)=ff(x(i,1),x(i,2));
 D(i)=F(i-1)-F(i);
 end
 if abs(x(n,:)-x0)<1e-4%收敛条件
 break;
 end
 d=x(n,:)-x0;
 x(n+1,:)=2*x(n,:)-x0;
 [D_m,j]=max(D);
 F2=F(n);F3=ff(x(n+1,1),x(n+1,2));
 if F3<F0 && (F0-2*F2+F3)*(F0-F2-D_m)^2>=D_m/2*(0-F3)^2
 x01=x(n,1);x02=x(n,2);
 opt_step=goldsection(x01,x02,d,h0); 
 opt_step;
 x(n+1,:)=x(n,:)+opt_step*d;
 e(j,:)=[];
 e(n,:)=d;
 x0=x(n+1,:);
 x_d(k,:) = x0; %%记录迭代路径
 else 
 if F2<F3
 x0=x(n,:);
 else
 x0=x(n+1,:);
 x_d(k,:) = x0; %%记录迭代路径
 end 
 end
end
disp('最优化结果 x0')
disp(x0)
F_x0 = ff(x0(1),x0(2));
disp('最优化结果最小函数值 F_x0')
disp(F_x0)
disp(k)
x_road = x_d(:,1);
y_road = x_d(:,2);
z_road = (1 - x_road).^2 + 100 .* (y_road - x_road.^2).^2;
xx = -0.5:0.05:1.2;
yy = -0.5:0.05:1.2;
[X,Y] = meshgrid(xx,yy);
Z = (1-X).^2 + 100 .* (Y - X.^2).^2 ;
figure
mesh(X,Y,Z);
xlabel('x')
ylabel('y')
zlabel('z')
hold on
plot3(x_road,y_road,z_road,'r*')

 子程序goldsection.m

function opt_step=goldsection(x01,x02,d,h0)
l=0.618; 
[a,b]=search2(x01,x02,d,h0); 
a1=b-l*(b-a);y1=ff(x01+d(1)*a1,x02+d(2)*a1);
a2=a+l*(b-a);y2=ff(x01+d(1)*a2,x02+d(2)*a2);
for n=1:100
 if y1>=y2
 x1(n)=a1;x2(n)=a2;yp1(n)=y1;yp2(n)=y2;
 a=a1;a1=a2;y1=y2;
 a2=a+l*(b-a);y2=ff(x01+d(1)*a2,x02+d(2)*a2);
 else
 x1(n)=a1;x2(n)=a2;yp1(n)=y1;yp2(n)=y2;
 b=a2;a2=a1;y2=y1;
 a1=b-l*(b-a);y1=ff(x01+d(1)*a1,x02+d(2)*a1);
 end
 
 aa(n)=(a+b)/2; 
 y(n)=ff(x01+d(1)*aa(n),x02+d(2)*aa(n));
 e(n)=abs(b-a);
 opt_step=(a+b)/2;
 if abs(b-a)<1e-5
 break;
 end
end

子程序 search2.m

function [a,b]=search2(x01,x02,d,h0) 
a1=0; y1=ff(x01+d(1)*a1,x02+d(2)*a1);
h=h0;
a2=h; y2=ff(x01+d(1)*a2,x02+d(2)*a2);
if y2>y1 
 h=-h;
 a3=a1;y3=y1;
 a1=a2;y1=y2;
 a2=a3;y2=y3;
end
a3=a2+h;y3=ff(x01+d(1)*a3,x02+d(2)*a3);
while y3<y2
 h=2*h;
 a1=a2;y1=y2;
 a2=a3;y2=y3;
 a3=a2+h; y3=ff(x01+d(1)*a3,x02+d(2)*a3); 
end
if h>0
 a=a1;b=a3;
else
 a=a3;b=a1;
end

目标函数ff.m

function y=ff(x1,x2) 
y =(1-x1)^2+100*(x2-x1^2)^2;

对于不同的初始点进行优化仿真得出结果如下：

 搜索路径如下图所示：