matlab的TPO算法怎么用

%Tree Physiology Optimization (TPO) 無約束尋優例程

function[Solution]=TPO_UnconstrainedOptim(n, Num_iterations)

tic;

if nargin<3, Num_iterations=30; end

if nargin<2, pop=30; end

if nargin<1, n=30; end

%目標函數開始（一個神奇的函數）

%--------------Shekel foxholes function---------------

strf1= '-((((x-9.681)^2+0.806)+ ((y-0.667)^2+0.806))^-1 + (((x-9.400)^2+0.517)+ ((y-2.041)^2+0.517))^-1 + (((x-8.025)^2+0.100)+ ((y-9.152)^2+0.100))^-1 + (((x-2.196)^2+0.908)+ ((y-0.415)^2+0.908))^-1 ';

strf2= '+(((x-8.074)^2+0.965)+ ((y-8.777)^2+0.965))^-1 + (((x-7.650)^2+0.669)+ ((y-5.658)^2+0.669))^-1 + (((x-1.256)^2+0.524)+ ((y-3.605)^2+0.524))^-1 + (((x-8.314)^2+0.902)+ ((y-2.261)^2+0.902))^-1 ';

strf3= '+(((x-0.226)^2+0.531)+ ((y-8.858)^2+0.531))^-1 + (((x-7.305)^2+0.876)+ ((y-2.228)^2+0.876))^-1 + (((x-0.652)^2+0.462)+ ((y-7.027)^2+0.462))^-1 +(((x-2.699)^2+0.491)+ ((y-3.516)^2+0.491))^-1 ';

strf4= '+(((x-8.327)^2+0.463)+ ((y-3.897)^2+0.463))^-1 + (((x-2.132)^2+0.714)+ ((y-7.006)^2+0.714))^-1 + (((x-4.707)^2+0.352)+ ((y-5.579)^2+0.352))^-1 + (((x-8.304)^2+0.869)+ ((y-7.559)^2+0.869))^-1'; 

strf5= '+(((x-8.632)^2+0.813)+ ((y-4.409)^2+0.813))^-1 + (((x-4.887)^2+0.811)+ ((y-9.112)^2+0.811))^-1 +(((x-2.440)^2+0.828)+ ((y-6.686)^2+0.828))^-1 + (((x-6.306)^2+0.964)+ ((y-8.583)^2+0.964))^-1'; 

strf6= '+(((x-0.652)^2+0.789)+ ((y-2.343)^2+0.789))^-1 + (((x-5.558)^2+0.360)+ ((y-1.272)^2+0.360))^-1 + (((x-3.352)^2+0.369)+ ((y-7.549)^2+0.369))^-1 + (((x-8.798)^2+0.992)+ ((y-0.880)^2+0.992))^-1'; 

strf7= '+(((x-1.460)^2+0.332)+ ((y-8.057)^2+0.332))^-1 + (((x-0.432)^2+0.817)+ ((y-8.645)^2+0.817))^-1 +(((x-0.679)^2+0.632)+ ((y-2.800)^2+0.632))^-1 + (((x-4.263)^2+0.883)+ ((y-1.074)^2+0.883))^-1'; 

strf8= '+ (((x-9.496)^2+0.608)+ ((y-4.830)^2+0.608))^-1 + (((x-4.138)^2+0.326)+ ((y-2.562)^2+0.326))^-1)';

funstr=strcat(strf1,strf2,strf3,strf4,strf5,strf6,strf7,strf8);

%目標函數終于結束了

% 把函數表達式字符串變成內聯函數再矢量化

f=vectorize(inline(funstr));

range= [-10 10 -10 10]; 

%劃分網格

NumOfGrid=50;

deltgridx=(range(2)-range(1))/NumOfGrid;

deltgridy=(range(4)-range(3))/NumOfGrid;

x_grid=range(1):deltgridx:range(2); y_grid=range(3):deltgridy:range(4);

[x,y]= meshgrid(x_grid,y_grid);

z=f(x,y);

surfc(x,y,z); 

%開始TPO算法

%初始化

[xn, yn, rnx, rny]=init_tree(pop,n);

figure(); grid on;                            

for jj=1:pop

    zn(jj,:)=f(xn(jj,:),yn(jj,:));          %初始目標函數值

    [z0(jj,1),Imin]=min(zn(jj,:));          %初始最優值

    x0(jj,1)=xn(jj,Imin);                  

    y0(jj,1)=yn(jj,Imin);                   

    r0x=rnx(jj,Imin);                      

    r0y=rny(jj,Imin);                     

end

[xn, yn, rnx, rny]=init_tree(pop,n);

for jj=1:n

    x00(:,jj)=x0;                          %初始枝干最優

    y00(:,jj)=y0;                         

end

z00=min(z0);                               %初始枝干最優值

xn=xn+(x00-xn)+cos(pi/3).*5.*rand(pop,n);  %植物生長

yn=yn+(y00-yn)+cos(pi/3).*5.*rand(pop,n);  

[xn, yn,range]=range_tree(xn, yn, range,pop);

%主要迭代過程

countt=1;                                   %記錄計數

for i=1:Num_iterations,  

    contour(x,y,z,15);                      %等高線圖

    hold on;                                

    %生長因數

    %alpha = 樹根伸長的吸收常數

    %beta = 樹枝的（經驗表明30和70之間較為合適）. 

    %theta = 降低隨機化的因素，隨著迭代增加而減少隨機化的冪數

    alpha=0.3;                                    

    beta= 50;            

    theta=0.9;           

    Theta=theta^(1*i);  

    %-----------------------------------------------------

    zn=f(xn,yn);

    for ij=1:pop

        [zb(ij,1),Imin]=min(zn(ij,:)); %當前最優

        xb(ij,1)= xn(ij,Imin);        

        yb(ij,1)=yn(ij,Imin);          

        if zb(ij,1)<z0(ij)             %和全局最優比較

            z0(ij)=zb(ij,1);

            x0(ij)=xb(ij,1);

            y0(ij)=yb(ij,1); 

        end

    end

    [z00,I00]=min(z0);                 %當前枝干最優

    for jj=1:n

        x00(:,jj)=x0(I00);             %當前整棵樹最優

        x0l(:,jj)=x0;                  %每一枝

        y00(:,jj)=y0(I00);             

        y0l(:,jj)=y0;                  

    end

    Cx=Theta*(x0l-xn);                 %生成養分

    Cy=Theta*(y0l-yn); 

    r0x=rnx; r0y=rny;                  

    rnx=rnx+ (rand(pop,n)-0.5).*alpha.*Cx;  %樹根生長

    rny=rny+(rand(pop,n)-0.5).*alpha.*Cy;

    Nx= Theta.* (rnx - r0x);                %吸收養分

    Ny= Theta.* (rny - r0y);                

    for ii=1:pop

        xn=xn+(x00-xn)+(Nx).*cos(ii*pi/3).* beta;  %樹枝生長

        yn=yn+(y00-yn)+(Ny).*cos(ii*pi/3).* beta; 

    end

    [xn, yn,range]=range_tree(xn, yn, range,pop);

    [x00, y00,range]=range_tree(x00, y00, range,pop);

    %每個樹枝達到最佳值
  

    plot(xn, yn,'.',x0, y0, '*','MarkerSize',5); axis(range);  

    drawnow;                                                   

    hold off;                                                     

    Solution(countt,:)= [z00(1,1)];

    countt=countt+1;

end

%初始化函數

function[xn, yn, rnx, rny]= init_tree(pop,n) %initialization value

xn=rand(pop,n);

yn=rand(pop,n);

rnx=rand(pop,n); 

rny=rand(pop,n);

%檢測邊界

function[xn, yn, range] = range_tree(xn, yn, range,pop)

nn=length(xn);

for ij=1:nn

    for iij=1:pop

        if xn(iij,ij)<=range(1), xn(iij,ij)=range(1); end

        if xn(iij,ij)>=range(2), xn(iij,ij)=range(2); end

        if yn(iij,ij)<=range(3), yn(iij,ij)=range(3); end

        if yn(iij,ij)>=range(4), yn(iij,ij)=range(4); end

    end