華電北風吹
天津大學認知計算與應用重點實驗室
日期：2015/6/28

粒子群演算法求解函式極值
初始化一個粒子群體隨機分佈在解空間中。然後每一次迭代，根據每個粒子的歷史最優位置和粒子群的最優位置更新每個粒子的位置，直到滿足要求為止。

clc;clear all;
tic;                              %程式執行計時
E0=0.001;                        %允許誤差
MaxNum=100;                    %粒子最大迭代次數
narvs=1;                         %目標函式的自變數個數
particlesize=30
 
;                    %粒子群規模
c1=2;                            %每個粒子的個體學習因子，也稱為加速常數
c2=2;                            %每個粒子的社會學習因子，也稱為加速常數
w=0.6;                           %慣性因子
vmax=0.8;                        %粒子的最大飛翔速度
x=-5+10*rand(particlesize,narvs);     %粒子所在的位置
v=2*rand(particlesize,narvs);         %粒子的飛翔速度
 

%用inline定義適應度函式以便將子函式檔案與主程式檔案放在一起，
%目標函式是：y=1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2))
%inline命令定義適應度函式如下：
fitness=inline('1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');
%inline定義的適應度函式會使程式執行速度大大降低
for i=1:particlesize
    for j=1:narvs
        f(i)=fitness(x(i,j));
    end
end
personalbest_x=x;
personalbest_faval=f;
[globalbest_faval i]
 
=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
k=1;
while k<=MaxNum
    for i=1:particlesize
        for j=1:narvs
            f(i)=fitness(x(i,j));
        end
        if f(i)<personalbest_faval(i) %判斷當前位置是否是歷史上最佳位置
            personalbest_faval(i)=f(i);
            personalbest_x(i,:)=x(i,:);
        end
    end
    [globalbest_faval i]=min(personalbest_faval);
    globalbest_x=personalbest_x(i,:);
    for i=1:particlesize %更新粒子群裡每個個體的最新位置
        v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...
            +c2*rand*(globalbest_x-x(i,:));
        for j=1:narvs    %判斷粒子的飛翔速度是否超過了最大飛翔速度
            if v(i,j)>vmax;
                v(i,j)=vmax;
            elseif v(i,j)<-vmax;
                v(i,j)=-vmax;
            end
        end
        x(i,:)=x(i,:)+v(i,:);
    end
    if abs(globalbest_faval)<E0,break,end
    k=k+1;
end
Value1=1/globalbest_faval-1; Value1=num2str(Value1);
% strcat指令可以實現字元的組合輸出
disp(strcat('the maximum value','=',Value1));
%輸出最大值所在的橫座標位置
Value2=globalbest_x; Value2=num2str(Value2);
disp(strcat('the corresponding coordinate','=',Value2));
x=-5:0.01:5;
y=2.1*(1-x+2*x.^2).*exp(-x.^2/2);
plot(x,y,'m-','linewidth',3);
hold on;
plot(globalbest_x,1/globalbest_faval-1,'kp','linewidth',4);
legend('目標函式','搜尋到的最大值');xlabel('x');ylabel('y');grid on;toc;

程式執行結果如下(參考適應度函式)