【MATLAB】POWELL法（鲍维尔法）求解二维优化问题

function powell
syms y x1 x2;
y=4*(x1-5)^2+(x2-6)^2;              %目标函数
x0=[8,9];
F0=0.00001;
F1=0;
F2=1;
F3=0;
fanwei=1:0.01:10;
[X,Y] = meshgrid(fanwei,fanwei);
Z=4*(X-5).^2+(Y-6).^2; 
contour(X,Y,Z,20);
hold on;
 while      abs(F2-F0)/F0>0.01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%进行第一个方向的搜索
plot(x0(1),x0(2),'ro');
x_temp=x0;
F0=subs(y,[x1,x2],x0);
step=1;
n=0;
x=[0 0 0];
% interval=[-100:0.01:150];
% plot(interval,f(interval))
% hold on
 flag=-1;   % 终止搜索标志
while subs(y,[x1,x2],[x0(1),x(1)])<=subs(y,[x1,x2],[x0(1),x(2)])|subs(y,[x1,x2],[x0(1),x(2)])>=subs(y,[x1,x2],[x0(1),x(3)])
    x(1)=x(2);
     x(2)=x(3);
    x(3)= x(3)+step*2^n;
    n=n+1;
   if abs(x(3))>10^10

   x=[0 0 0];
   step=-1;
   n=0;
  
  flag=flag+1;
   if flag==1
       
       break;
   end
   end
  
end
%  plot(x,f(x),'o');
 a=x(1);b=x(3);
x1temp=0.382*(b-a)+a;
x2temp=0.618*(b-a)+a;
while (b-a)/b>0.0002
if subs(y,[x1,x2],[x0(1),x1temp])<subs(y,[x1,x2],[x0(1),x2temp])
    b=x2temp;
    x2temp=x1temp;
    x1temp=0.382*(b-a)+a;
else
    a=x1temp;
    x1temp=x2temp;
    x2temp=0.681*(b-a)+a;
    
end
end
if flag~=1
xo=(a+b)/2 ; 
x_temp(2)=xo;
F1=subs(y,[x1,x2],x_temp);
end

%%%%%%%%%%%进行第二个方向的搜索
 plot(x_temp(1),x_temp(2),'ro');
step=.1;
n=0;
x=[0 0 0];
% interval=[-100:0.01:150];
% plot(interval,f(interval))
% hold on
 flag=-1;
while subs(y,[x1,x2],[x(1),xo])<=subs(y,[x1,x2],[x(2),xo])|subs(y,[x1,x2],[x(2),xo])>=subs(y,[x1,x2],[x(3),xo])
    x(1)=x(2);
     x(2)=x(3);
    x(3)= x(3)+step*2^n;
    n=n+1;
   if abs(x(3))>10^5
   x=[0 0 0];
   step=-1;
   n=0;
     
  flag=flag+1;
   if flag==1
       break;
   end
   end
  
end
%  plot(x,f(x),'o');
 a=x(1);b=x(3);
x1temp=0.382*(b-a)+a;
x2temp=0.618*(b-a)+a;
while (b-a)/b>0.0002
if subs(y,[x1,x2],[x1temp,xo])<subs(y,[x1,x2],[x2temp,xo])
    b=x2temp;
    x2temp=x1temp;
    x1temp=0.382*(b-a)+a;
else
    a=x1temp;
    x1temp=x2temp;
    x2temp=0.681*(b-a)+a;
    
end
end
if flag~=1
xo=(a+b)/2;
x_temp(1)=xo;
end

F2=subs(y,[x1,x2],x_temp);
x_flec=2*x_temp-x0;              %反射点

F3=subs(y,[x1,x2],x_flec);
deta1=F0-F1;
deta2=F1-F2;
detam=max([deta1,deta2]);
if F3<F0  &  (F0-2*F2+F3)*(F0-F2-detam)^2<0.5*detam*(F0-F3)^2
d=x_temp-x0;
rate=d(2)/d(1);


%进行d方向求最优点
step=1;
n=0;
x=[0 0 0];
% interval=[-100:0.01:150];
% plot(interval,f(interval))
% hold on
 flag=-1;
while subs(y,[x1,x2],[x(1),rate*x(1)])<=subs(y,[x1,x2],[x(2),rate*x(2)])|subs(y,[x1,x2],[x(2),rate*x(2)])>=subs(y,[x1,x2],[x(3),rate*x(3)])
    x(1)=x(2);
     x(2)=x(3);
    x(3)= x(3)+step*2^n;
    n=n+1;
   if abs(x(3))>10^5
   x=[0 0 0];
   step=-1;
   n=0;
  
  flag=flag+1;
   if flag==1
       break;
   end
   end
  
end
%  plot(x,f(x),'o');
 a=x(1);b=x(3);
x1temp=0.382*(b-a)+a;
x2temp=0.618*(b-a)+a;
while (b-a)/b>0.0002
if subs(y,[x1,x2],[x1temp,rate*x1temp])<subs(y,[x1,x2],[x2temp,rate*x2temp])
    b=x2temp;
    x2temp=x1temp;
    x1temp=0.382*(b-a)+a;
else
    a=x1temp;
    x1temp=x2temp;
    x2temp=0.681*(b-a)+a;
    
end
end
if flag~=1
xo=(a+b)/2;
x0=[xo,xo*rate];
end

else
if F3<F2
    x0=x_flec;
else
    x0=x_temp;
end
end
end

x0
plot(x0(1),x0(2),'r*');
end