github: 智慧演算法的課件和參考資料以及實驗程式碼

粒子群演算法的原理:

粒子群演算法是一種群體智慧演算法，通過追隨當前搜尋到的最優值來尋找全域性最優。該演算法實現容易、精度高、收斂快，在解決實際問題中具有很大的優越性。主要步驟可描述如下：

1、初始化粒子群位置和速度。

2、計算每個粒子的適應度，確定全域性最優粒子gbest和個體最優粒子pbest。

3、判斷演算法收斂準則是否滿足，若滿足，則輸出搜尋結果；否則執行以下步驟。

4、根據gbest和pbest更新速度分量。

5、根據更新的速度分量更新粒子位置。

6、更新全域性最優粒子和個體最優粒子。

7、返回步驟③。

下面我們舉了兩個函式最小值優化問題，實際上是引數擬合，但是轉化成了最小化問題(具體下載第三次實驗課件)

f.m

function y = f( x )

% 矩陣按照元素的平方是X.^2 A*B 是矩陣的乘法，而A.*B是兩個矩陣對應元素相乘，其他情況下效果一樣
% X=[1,3,4,5,6,7,8,9,10]; 
% Y=[10,5,4,2,1,1,2,3,4];
% Y1=x(1)*X.^2+x(2)*X+x(3);

t=1:15;
Y=[352	211	197	160	142	106	104	60	56	38	36	32	21	19	15];
Y1=x(1)*exp(x(2).*t);

y=norm(Y-Y1); % 向量的歐幾里得距離

end

 

gbestpso.m

clear
clc
%*********初始化
M=20;  %種群規模
x1=rand(M,1); % 初始化粒子位置
x2=rand(M,1);
%x3=rand(M,1);
%X=[x1,x2,x3]
X=[x1,x2];
c1=2; % c1和c2是學習因子
c2=2;
wmax=0.9; % 最大最小慣性權重
wmin=0.4; % 線性減小慣性權重
Tmax=50; % 迭代次數
%v=zeros(M,3);
v=zeros(M,2); % 初始化速度
% v1m=0.1;  % 速度1約束
% v2m=1;% 速度2約束
%******* 全域性最優粒子位置初始化
fmin=1000; 
for i=1:M
    fx = f(X(i,:));
    if fx<fmin
        fmin=fx;
        gb=X(i,:);
    end
end
%********粒子個體歷史最優位置初始化
pb=X; 
%********演算法迭代
for t=1:Tmax
    t % w是每代對應的慣性權重，正在逐漸減小
    w=wmax-(wmax-wmin)*t/Tmax;  % 線性下降慣性權重
    for i=1:M
       %******更新每個粒子的當前速度
       v(i,:)=w*v(i,:)+c1*rand(1)*(pb(i,:)-X(i,:))+c2*rand(1)*(gb-X(i,:));
       %******速度約束
%        if v(i,1)>v1m; v(i,1)=v1m;end
%        if v(i,2)>v2m; v(i,2)=v2m; end  
%        if v(i,1)<-v1m;v(i,1)=-v1m;end
%        if v(i,2)<-v2m;v(i,2)=-v1m;end
       %*******更新每個粒子的當前位置
       X(i,:)=X(i,:)+v(i,:);
    end
    % 更新所有粒子的速度和位置之後更新pbest和gbest
    for i=1:M
        if f(X(i,:))<f(pb(i,:))
            pb(i,:)=X(i,:);
        end
        if f(X(i,:))<f(gb)
            gb=X(i,:);
        end
    end
    % 儲存最佳適應度
    re(t)=f(gb);
end
figure(1)
plot(re)
xlabel('迭代次數')
ylabel('適應度函式')
gb
re(end)
figure(2)
% X=[1,3,4,5,6,7,8,9,10]; 
% Y=[10,5,4,2,1,1,2,3,4];
% Y1=gb(1)*X.^2+gb(2)*X+gb(3);
X=1:15;
Y=[352  211 197 160 142 106 104 60  56  38  36  32  21  19  15];
Y1=gb(1)*exp(gb(2).*X);

t=16;
y16=gb(1)*exp(gb(2).*t)
plot(X,Y,'o',X,Y1,':');
hold on
plot(t,y16,'ro')

% gb = [400.2029 -0.2241] 最好適應度61.2758 
 

lbestpso.m

clear
clc
%*********初始化
M=20;  % 種群規模
x1=rand(M,1); % 初始化例子位置
x2=rand(M,1);
%x3=rand(M,1);
%X=[x1,x2,x3]
X=[x1,x2];
c1=2;  % 學習因子
c2=2;
wmax=0.9; % 最大最小慣性權重
wmin=0.4;
Tmax=50; % 迭代次數
v=zeros(M,2); % 初始化速度
%*******全域性最優粒子位置初始化
fmin=1000;
gb=X;
k=2; % 鄰域規模
for i=1:M
    for in=i-k:i+k
        if in<1;
            in=in+M;%越界處理
        elseif in>M
            in=in-M;%越界處理
        else
            in=in;
        end
        if f(X(in,:))<f(gb(i,:))
            gb(i,:)=X(in,:);%個體的鄰域最優
        end
    end
end
%********粒子個體歷史最優位置初始化
pb=X; 
%********演算法迭代
for t=1:Tmax
    t
    w=wmax-(wmax-wmin)*t/Tmax;  %線性下降慣性權重
    for i=1:M
       %******更新粒子速度 終於是區域性最優
       v(i,:)=w*v(i,:)+c1*rand(1)*(pb(i,:)-X(i,:))+c2*rand(1)*(gb(i,:)-X(i,:));
       %*******更新粒子位置
       X(i,:)=X(i,:)+v(i,:);
    end
    %更新pbest和gbest
    for i=1:M
        if f(X(i,:))<f(pb(i,:))
            pb(i,:)=X(i,:);
        end
        for in=i-k:i+k
            if in<1;
                in=in+M;
            elseif in>M
                in=in-M;
            else
                in=in;
            end
            if f(X(in,:))<f(gb(i,:))
                gb(i,:)=X(in,:);%更新鄰域最優
            end
        end
    end
    %儲存最佳適應度
    ffmin=1000;
    for i=1:M
        if f(gb(i,:))<ffmin;
            ffmin=f(gb(i,:))
            re(t)=ffmin;
            ggb=gb(i,:);
        end
    end
end
figure(1)
plot(re)
xlabel('迭代次數')
ylabel('適應度函式')
ggb
figure(2)
%X=[1,3,4,5,6,7,8,9,10]; 
%Y=[10,5,4,2,1,1,2,3,4];
%Y1=gb(1)*X.^2+gb(2)*X+gb(3);
X=1:15;
Y=[352 211  197 160 142 106 104 60  56  38  36  32  21  19  15];
Y1=ggb(1)*exp(ggb(2).*X);
t=16;
y16=ggb(1)*exp(ggb(2).*t)
plot(X,Y,'o',X,Y1,':');
hold on
plot(t,y16,'ro')
hold off
 
% 最好適應度61.2860 指數函式引數a,b分別是400.0875 -0.2234 
% 二次函式引數擬合問題 a = -0.0963 b = 0.7585 c = 3.2129 最小適應度函式為8.5722

下面是在matlab上執行的截圖:

lbest二元函式引數擬合圖:

lbest指數函式引數擬合:

gbest指數函式引數問題: