基于黄金分割的修正Powell算法

1.目标函数

function y= TextF(x)
%目标函数
%输入参数x[x1,x2];
%输出参数y
y = x(1)*x(1)+2*x(2)*x(2)-4*x(1)-2*x(1)*x(2);
end</span>

2.修订的Powell算法

function [Matrix,Value] = Powell
%--------------------------------Step1-------------------------------------
%设置初始迭代轮次k=1；迭代精度：Err；初始点；初始搜索方向组；
Err = 0.01;
X0 = [1,1];
Dir = [1,0;
       0,1];
num = 2; 

while(1)
%----------------------------------Step2-----------------------------------
%第一轮：
    %沿着第一个方向进行一维搜索的得到相应的极值
    [X(1,:),JZ(1)] = OneDimensionSearch(X0,Dir(1,:),Err);
    %沿着第2到第num方向进行一维搜索
    for k = 2:num
        dir(k,:) = Dir(k,:);
        [X(k,:),JZ(k)] = OneDimensionSearch(X(k-1,:),dir(k,:),Err);
    end
    %构造新的搜索方向向量 并以X0为起点，沿此方向进行新一轮搜索，得到相应的极值
    dir(num+1,:) = X(num,:)-X0;
    [X(num+1,:),JZ(num+1)] = OneDimensionSearch(X0,dir(num+1,:),Err);
%----------------------------------Step3-----------------------------------
    %判断迭代终止条件
    delta = X(num+1,:)-X0;
    if( sqrt(sum(delta.*delta))<=Err )
        Matrix = X(num+1,:);
        Value = TextF(Matrix); %输出最优解和最优变换参数
        break;
    else
%----------------------------------Step4-----------------------------------
        MuInfoVal_0 = TextF(X0);
        %计算各个方向上的下降量
        Diff(1) = MuInfoVal_0-JZ(1);%此时我们在寻找极小值
        for k = 2:num
            Diff(k) =  JZ(k-1)-JZ(k);
        end
        %记录最下降量 以及 其出现的位置
        [maxDiff,m] = max(Diff);
%----------------------------------Step5----------------------------------- 
       %计算映射点
       Xmap = 2*X(num,:)-X0;
       %计算Powell特征值
       Val_0 = TextF(X0);
       Val_1 = TextF(X(num,:));
       Val_2 = TextF(Xmap);
       %判断Powell条件
       if(Val_2>Val_0 && ...
          (Val_0-2*Val_1+Val_2)*(Val_0-Val_1-maxDiff)<0.5*maxDiff(Val_0-Val_2)^2 )
%----------------------------------Step6----------------------------------- 
       %如果满足,即为非线性相关 %设置下一轮搜索的初始点
       X0 = X(num+1,:);
       %替换m位置的搜索方向
       Dir(m,:) = dir(num+1,:);
%----------------------------------Step7----------------------------------- 
       else %如果不满足，即为线性相关，%保持原来的搜索方向，设置下一轮搜索始点
           if(Val_1 < Val_2)
               X0 = X(num,:);
           else
               X0 = X(num+1,:);
           end
       end 
    end    
end
end

3.一维搜索策略（黄金分割、数值解法）

%% 一维搜索函数:采用黄金分割法 %%
function [ExtremePos,ExtremeVal] = OneDimensionSearch(x1,Dir,err)
% 从X0点开始,进行搜索
% Dir:搜索的方向。
% err:黄金分割的停止搜索的精度。
% 返回极大值ExtremeVal，以及极大值对应的点ExtremePos
%------------------------Step1：进退法确定区间------------------------------
%计算x1,前进，计算x2
y1 = TextF(x1);
x2 = x1+Dir;
y2 = TextF(x2);
%最小值左右支判断
if y1 < y2 % 右支
    Dir = -Dir;
    temp = x1;
    x1 = x2;
    x2 = temp;
    x3 = x2+Dir;
    y3 = TextF(x3);
else %左支 扩大搜索步法
    Dir = 2*Dir;
    x3 = x2+Dir;
    y3 = TextF(x3);
end
while(1) 
    if (y2<=y3) %发现极小值所在的区间
        a = min(x1,x3);
        b = max(x1,x3);
        break;
    else
        x1 = x2;
        x2 = x3;
        y2 = y3;
        x3 = x2+Dir;
        y3 = TextF(x3);
    end
end
%--------------------------Step2：黄金分割----------------------------------
g1 = 0.382; g2 = 0.618;
xx1 = a+g1*(b-a);
yy1 = TextF(xx1);
xx2 = a+g2*(b-a);
yy2 = TextF(xx2);
delta = sqrt(sum(abs(b-a).*abs(b-a)));
while( delta>=err)
    if(yy1 < yy2)
        b = xx2;
        xx2 = xx1;
        yy2 = yy1;
        xx1 = a+g1*(b-a);
       yy1 = TextF(xx1);
    else
        a = xx1;
        xx1 = xx2;
        yy1 = yy2;
        xx2 = a+g2*(b-a);
        yy2 = TextF(xx2);
    end   
    delta = sqrt(sum(abs(b-a).*abs(b-a)));
end
ExtremePos = 0.5*(a+b);
ExtremeVal = TextF(ExtremePos);
end

4.运行代码及结果

[m,n]=Powell
m =
    3.9969    1.9985

n =
   -8.0000</span>