matlab鲍威尔法（Powell）

%改进鲍威尔法
clc
clear all

syms x1 x2 a1 a2 a3 b1 b2 b3
f = x1*x1+2*x2*x2-4*x1-2*x1*x2;
X=[x1;x2];
X0=[1 1]';   %给定初始点
Theta_error=0.001;
e1=[1 0]';
e2=[0 1]';

FX0=subs(f,X,X0);

n=2;
k=0;
k1=0;
k2=0;
d1=e1;%初始搜搜方向
d2=e2;
i=1;

%求a1，a2，X1,X2
Da1=[];
Da2=[];
Da3=[];
DX0=[];
DX00=[];
DXX0=[];
DXXX0=[];
DX1=[];
DX11=[];
DX2=[];
DX22=[];
DX3=[];
DX33=[];
DF=[];

Xh1 = X0+a1.*d1;
FXh1 =  subs(f,X,Xh1);
dFXh1 = gradient(FXh1,a1);
a1 = solve(dFXh1);
f = x1*x1+2*x2*x2-4*x1-2*x1*x2;

X1 = X0+a1*e1;
FX1 = subs(f,X,X1);

Xh0 = X1+a2.*d2;
FXh2 = subs(f,X,Xh0);
dFXh0 = gradient(FXh2,a2);
a2 = solve(dFXh0);
X2 = X1+a2*e2;
FX2 = subs(f,X,X2);

deta1 = FX0-FX1;
deta2 = FX1-FX2;
m = [deta1,deta2];
M=sort(m);   %将deta从小到大排序，然后最后一位就是最大值
deta_Max=M(2);

%k=1
d3=X2-X0;
X3=2*X2-X0;
FX3 = subs(f,X,X3);

if FX3 < FX0 && (FX0-2*FX2+FX3)*(FX0-FX2-deta_Max)^2 < 0.5*deta_Max*(FX0-FX3)^2
    XX0=X2+a3.*d3;
    h0=X2+a3.*d3;
    Fh0 = subs(f,X,h0);
    dFh0 = gradient(Fh0,a3);
    a3 = solve(dFh0);    
    XX0=X2+a3.*d3;
    d1=d2;
    d2=d3;
else
    d1=d1;
    d2=d2;
    if FX2<FX3
        XX0=X2;
    else
        XX0=X3;
    end
end
while  norm(XX0-X0) >= Theta_error
    
    k=k+1;
    X0=XX0;  DX0(k)=X0(1);DX00(k)=X0(2);
    %求a1，a2，X1,X2
    a1=b1;
    Xh1 = XX0+a1.*d1;
    FXh1 =  subs(f,X,Xh1);
    dFXh1 = gradient(FXh1,b1);
    a1 = solve(dFXh1);    Da1(k)=a1;
    f = x1*x1+2*x2*x2-4*x1-2*x1*x2;
    
    
    X1 = XX0+a1.*d1;   DX1(k)=X1(1);DX11(k)=X1(2);
    FX1 = subs(f,X,X1);
    
    a2=b2;
    Xh0 = X1+a2.*d2;
    FXh2 = subs(f,X,Xh0);
    dFXh0 = gradient(FXh2,b2);
    a2 = solve(dFXh0);  Da2(k)=a2;
    X2 = X1+a2.*d2;  DX2(k)=X2(1);DX22(k)=X2(2);
    FX2 = subs(f,X,X2);
    
    FX0 = subs(f,X,XX0);
    
    deta1 = FX0-FX1;
    deta2 = FX1-FX2;
    m = [deta1,deta2];
    M=sort(m);   %将deta从小到大排序，然后最后一位就是最大值
    deta_Max=M(2);
    
    
    if FX3 < FX0 && (FX0-2*FX2+FX3)*(FX0-FX2-deta_Max)^2 < 0.5*deta_Max*(FX0-FX3)^2
        %走该条通道看k1++
        k1=k1+1;
        a3=b3; Da3(k)=a3;
        XX0=X2+a3*d3;
        h0=X2+a3*d3;
        Fh0 = subs(f,X,h0);
        dFh0 = gradient(Fh0,a3);
        a3 = solve(dFh0);  Da3(k)=a3;
        XX0=X2+a3*d3; DXX0(k)=XX0(1);DXXX0(k)=XX0(2);DF(k)=subs(f,X,XX0);
        d1=d2;
        d2=d3;
    else
              %走该条通道看k2++
        k2=k2+1;
        if FX2<FX3
            XX0=X2; DXX0(k)=XX0(1);DXXX0(k)=XX0(2);DF(k)=subs(f,X,XX0);
        else
            XX0=X3; DXX0(k)=XX0(1); DXXX0(k)=XX0(2);DF(k)=subs(f,X,XX0);
        end
    end
    
    %axis([0 10 0 10]);
    plot(XX0(1),XX0(2),'rp')
    %plot(X0(1),X0(2),'b*')
    line([X0(1) X1(1)],[X0(2) X1(2)]);
    line([X0(1) X2(1)],[X0(2) X2(2)]);
    line([X1(1) X2(1)],[X1(2) X2(2)]);
    line([X0(1) XX0(1)],[X0(2) XX0(2)]);
    hold on
end
grid on
axis([3 5 0 4])
ezplot('-8 =x1*x1+2*x2*x2-4*x1-2*x1*x2')
hold on
ezplot('-7.998 =x1*x1+2*x2*x2-4*x1-2*x1*x2')
hold on
ezplot('0 =x1*x1+2*x2*x2-4*x1-2*x1*x2')
hold on

X_output=XX0;
FX_output= subs(f,X,X_output);

%double(X_output);double(FX_output);
X_output,FX_output