黄金分割に基づく修正Powellアルゴリズム

1.目標関数

<span style="font-size:14px;">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アルゴリズム

f<span style="font-size:14px;">unction [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</span>

3.一次元検索ポリシー(金分割、数値解法)

<span style="font-size:14px;">%%       :        %%
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</span>

4.運転コード及び結果

<span style="font-size:14px;">[m,n]=Powell
m =
    3.9969    1.9985

n =
   -8.0000</span>