為了求解優化問題，有很多生物啟發算法。通過生物的某些特性來避免局部最小、加快尋優。這里要介紹的是香港理工大學（The Hong Kong Polytechnic University）的Yinyan Zhang 和 Shuai Li 提出的Porcellio scaber算法（PSA），學習的應該是中學生物書上講的鼠婦的生存規則。

function PS_simple(number_of_PS,MaxStep)

% porcellio scaber的數量

N=number_of_PS; 

%initalize a matrix to store position data 

x=zeros(MaxStep,2,N); 

%Generate the initial positions of all the porcellio scaber

x(1,:,:)=4*rand(2,N); 

%Set weighted paramter \lambda for decion based on aggregation and the

%propensity to explore novel enviroments

lambda=0.8;

%generate a series of \tau and make each elment of \tau a zero mean

%random real number 

sigma=0.001; %standard deviation 

tau_data=zeros(MaxStep,2,N);

for i=1:2

    for j=1:N

    tau_data(:,i,j)=sigma*(2*rand(MaxStep,1)-1);

    end

end

iter=1;

while iter<MaxStep

    %Get the position with the best environment condition at the current

    %time among the group of porcellio scaber

    minf=min(fun(x(iter,1,:),x(iter,2,:)));

    f=fun(x(iter,1,:),x(iter,2,:));

    [~,indmin]=find(abs(f-minf)<eps);

    % x-axis coordinate of the current best position 

    x_best_x=x(iter,1,indmin); 

    % y-axis coordinate of the current best position

    x_best_y=x(iter,2,indmin);  

    %best current position 

    x_best=[x_best_x,x_best_y]; 

    %Randomly choose a direction \tau to detect 

    all_tau=tau_data(iter,:,:);

    %detect the best enviroment condition minEx and the worst environment

    %condition maxEx at position \mathbf{x}^k_i+\tau for i=1:N all N

    %porcellio scaber

     x_p_tau=zeros(1,2,N); 

     for i=1:N

         %calculate \mathbf{x}^k_i+\tau

         x_p_tau(:,:,i)=x(iter,:,i)+all_tau(:,:,i);

     end

     Ex_p_tau=zeros(N,1);

     for i=1:N

         %calculate f(\mathbf{x}^k_i+\tau)

         Ex_p_tau(i)=fun(x_p_tau(:,1,i),x_p_tau(:,2,i)); 

     end

     %get max{f(\mathbf{x}^k_i+\tau)}

     maxEx_p_tau=max(Ex_p_tau); 

     %get min{f(\mathbf{x}^k_i+\tau)}

     minEx_p_tau=min(Ex_p_tau); 

    for i=1:N

        x_k_i=x(iter,:,i);

        %Determine the difference with respect to the position to aggregate

        diff=x_k_i-x_best;

        %Determine where to explore

        p_tau=calculate_p_tau(i,x_k_i,maxEx_p_tau,minEx_p_tau,all_tau);

        %movement according to the weighted result of aggregation and the 

        %propensity to explore novel enviroments

        x(iter+1,:,i)=x_k_i-(1-lambda)*diff-lambda*p_tau; 

    end

    iter=iter+1; %update iteration number

end

display('The optimal solution is: ')

x_best

display('The corresopnding function value is: ');

fun(x_best(1),x_best(2))

%visualization

pseudoFig=figure;

draw_pro; %draw the pseudo color figure of the problem

hold on;

%draw the tracjetory of all the procellio scaber

for i=1:N

     data_x=x(:,1,i);

     data_y=x(:,2,i);

plot(data_x,data_y,'k-.','LineWidth',1);

hold on

end

saveas(pseudoFig,'ex1result','fig');

%save all the data

save ALL_data 

save parameter_setting number_of_PS MaxStep

%The following are two subfunctions

%calculation of direction to explore

function p_tau=calculate_p_tau(i,x_k_i,maxEx_p_tau,minEx_p_tau,all_tau)

Ex_k_i_p_tau=(fun(x_k_i(1)+all_tau(1,1,i),x_k_i(2)+all_tau(1,2,i)));

p_tau=(Ex_k_i_p_tau-minEx_p_tau)/(maxEx_p_tau-minEx_p_tau)*all_tau(1,:,i);

end

end

通過在命令行執行PS_simple(20,40);運行

找到最小值后顯示

這是要求解的目標函數

function yout=fun(x,y)

yout=-sin(x).*(sin(x.^2/pi)).^20-sin(y).*(sin(2*y.^2/pi)).^20;

end

這里有一個將三維圖顯示在二維上的方法：以顏色深度表示本來的Z軸

[x,y]=meshgrid(linspace(0,4));

h=pcolor(x,y,fun(x,y));

set(h,'edgecolor','none','facecolor','interp');

colorbar;