粒子群最適化(PSO)とTCPSOのMatlabプログラム

粒子群最適化(PSO:Particle Swarm Optimization)

PSOは，Kennedy and Eberhartによって開発された，魚や鳥などにみられる群行動の餌の在り処の情報が群全体に伝達され，他の個体の行動に合わせて速度を調整しながらすべての個体が一斉に行動しているという点に着目して提案された最適解探索のアルゴリズムである．

PSOでは，位置情報と速度情報を持った粒子(個体)が最適解を求めて探索空間を探索する．各粒子は自分自身の過去の情報と周りの粒子の過去の情報を互いに共有しながら解を探索する．粒子の更新が単純な概念から容易に実装でき，さらに多様な改良が可能である．

粒子の更新

粒子$i,i=1,2,…,n$は時刻$t$における位置情報ベクトル$x_i(t)=(x_i1,x_i2,…,x_m)$と速度ベクトル$v_i (t)=(v_i1,v_i2,…,v_m)$の情報を持っている．

このとき，$n$は粒子数，$m$は次元数である．粒子$i$が各時刻まで探索した解のうち最良の位置情報ベクトルを$P_i(t)$として保持し，粒子群の中で最良の粒子の位置情報ベクトルを$P_G(t)$として群全体で共有する．

以下に位置情報ベクトルと速度ベクトルの更新式を示す．

$$ v_i(t+1)=w{v}_i(t)+C_1r_1(P_i(t)-m{x}_i(t))+C_2r_2(P_G(t)-x_i(t)) $$

$$ x_i(t+1)=x_i(t)+{v}_i(t+1) $$

$w$は慣性項に対する重み，$C_1$は自身の最良位置に対する重み，$C_2$は全体の最良位置に対する重みである．

また，$r_1,r_2∈[0,1]$は各次元ごとに毎期生成される一様分布に従ってランダムに決定される実数値である．

更新の概念図を以下に示す．

アルゴリズム

PSOのアルゴリズムを以下に示す．

粒子の初期生成初期の粒子の位置情報ベクトル$x_i$と速度ベクトル$v_i$を乱数に基づいて生成する．
評価値の算出各粒子の位置から評価値を決定する．
探索開始($t=0$)から現在までの各粒子$i$の最良位置情報$P_G (t)$の更新，各粒子ごとに，現在までの最良の評価値を記録する．
群全体の最良位置情報$P_G (t)$の更新全粒子の中から，最良の評価値を記憶する．
各粒子の速度を計算
各粒子の位置を計算
終了条件判定終了条件を満たすまで，Step2からStep6を繰り返す．

アルゴリズムの性質

単純な算術演算により構成され，パラメータ$w，C_1，C_2$への依存性が高いことや，各粒子が群全体の最良位置情報$P_G (t)$の位置に集中してしまうことで局所解に陥り，十分な探索ができないといった問題が存在する．収束速度と多様性の維持のバランスが難しい．

以下に，粒子数100個をランダムに生成し，$(x,y)=(0,0)$からの距離の逆数を評価値としてPSOを行った例を示す．

$t=0$のとき

$t=250$のとき

$t=500$のとき

収束は早いが，局所解に陥っていることがわかる．また，慣性項の影響で集中的に探索することが難しくなっている．

TCPSO(Two-swarm Cooperative Particle Swarm Optimization)

TCPSOでは大域的な探索を行うマスター粒子群と集中的な探索を行うスレーブ粒子群を用いることで，PSOの収束速度と多様性の問題を改善する．2つの粒子群がお互いに情報を共有しながら各粒子が位置と速度の更新を行う．性質の異なる粒子群を用いることで大域探索と集中探索を行うことが可能になり，効率的な解探索が可能になる．

粒子の更新

スレーブ粒子を添字$S$，マスター粒子を添字$M$で表す．以下に位置情報ベクトルと速度ベクトルの更新式を示す．

$$ v_i^S(t+1)=C_1^Sr_1(P_i^S(t)-x_i^S(t))+C_2^Sr_2(P_G(t)-x_i^S(t)) $$ $$ x_i^S(t+1)=x_i^S(t)+v_i^S(t+1) $$ $$ v_i^M(t+1)=w^Mv_i^M(t)C_1^Mr_1(P_i^M(t)-x_i^M(t)) +C_2^Mr_2(P_G^S(t)-x_i^M(t))+C_3^Mr_3(P_G(t)-x_i^M(t)) $$ $$ x_i^M(t+1)=x_i^M(t)+v_i^M(t+1) $$ マスター粒子の速度の更新には，スレーブ粒子群の全体の最良位置$P_G^S (t)$が追加されており，スレーブ群の情報を更新することで多様性を保持しながら大域的な探索をすることが可能になっている．スレーブ粒子とマスター粒子の更新の概念を以下に示す．

スレーブ粒子群の更新の概念図
マスター粒子群の更新の概念図

アルゴリズム

TCPSOのアルゴリズムを以下に示す．

粒子の初期生成スレーブ粒子群，マスター粒子群それぞれの初期の粒子の位置情報ベクトル $x_i^S,x_i^M$と速度ベクトル$v_i^S,v_i^M$を事前に決めた範囲よりランダムに生成する．
評価値の算出各スレーブ粒子と各マスター粒子の位置から評価値を決定する．
探索開始($t=0$)から現在までの各粒子の最良位置情報$P_i^S(t),P_i^M (t)$の更新，各粒子ごとに，現在までの最良の評価値を記憶する．
最良位置情報の更新スレーブ群全体での最良位置情報$P_G^S (t)$と，全体での最良位置情報$P_G (t)$の更新を記憶する．
各スレーブ群の粒子の速度を計算
各マスター群の粒子の速度，位置を計算
終了条件終了条件を満たしていれば終了し，満たさなければStep2に戻る．

アルゴリズムの性質

TCPSOはスレーブ群の収束性から，単峰型の問題に対しては基本的なPSOよりも早く収束させることができる．また，マスター群は多様性を保持しており，単峰型の問題の局所解への対応も可能となる．しかし，高次元，広範囲の問題に対しては，局所解に収束してしまうため多様性の維持が困難になり，適切な探索を行うことが難しい．

以下に，スレーブ粒子100個，マスター粒子100個をランダムに生成し，$(x,y)=(0,0)$からの距離の逆数を評価値としてPSOを行った例を示す．オレンジがスレーブ粒子群，ブルーがマスター粒子群である．

$t=0$のとき
$t=250$のとき
$t=500$のとき

スレーブ粒子群は，慣性項がないため収束が早くより集中的な探索ができるようになっていることがわかる．収束速度と多様性が確保されている．高次元で複雑非線形な問題に対しても，マスター粒子群が局所解を回避，逸脱し，多様な探索をし続けることが重要となる．

プログラム

PSO

%/////// <PSOパラメータ定義> /////////
ns=300;          % 粒子数
c1s=1;           % 自身の最良位置に対する重み
c2s=1;           % 全体の最良位置に対する重み
maxT=500;     % 繰り返し回数
D=2;             % 次元

%****************　[STEP1] 初期粒子群の発生 ********************************
            
ps = 0.5-rand(ns,2);  %初期粒子群発生 
ps(:,3) = 0;
ps_v = rand(ns,2);  
ps_best = ps; %ローカルベスト   
P_S_G(1,D+1) = 0; %グローバルベスト

f0 = figure;
plot(ps(:,1),ps(:,2),'o')
xlim([-0.5 0.5])
ylim([-0.5 0.5])

for count=1:maxT
%****************　[STEP2] 評価値計算 **************************************
            
    for t=1:ns
        ps(t,3) = 1/abs(ps(t,1))+abs(ps(t,2));
    end
               
%****************　[STEP3] 各粒子の最良位置の更新 ***************************
            
    for i=1:ns 
        if(ps(i,3) > ps_best(i,3))
            for j=1:D+1
               ps_best(i,j) = ps(i,j);
            end
        end
    end            
           
 %****************　[STEP4] 全体の最良位置の更新 ***************************
            
    for i=1:ns 
        if(ps(i,3) > P_S_G(1,3))
            for j=1:D+1
                P_S_G(i,j) = ps(i,j);
            end
        end 
    end
                        
 %****************　[STEP5] 各粒子の速度，位置を計算 ***************************     
            
    %パラメータ更新              
    for i=1:ns
        for j=1:D
            r1 = rand/100;
            r2 = rand/100;
            ps_v(i,j) = c1s*r1*(ps_best(i,j)-ps(i,j)) + c2s*r2*(P_S_G(1,j)-ps(i,j)); %速度更新
            ps(i,j) = ps(i,j) + ps_v(i,j); %位置更新
        end
    end

    if (mod(count,250)==0) %250更新ごとに図を表示
        f1 = figure;
        plot(ps(:,1),ps(:,2),'o')
        xlim([-0.5 0.5])
        ylim([-0.5 0.5])
    end
               
end

TCPSO

%/////// <TCPSOパラメータ定義> /////////
ns=300;          % スレーブ粒子数
nm=300;          % マスター粒子数
w=0.05;          %　慣性項
c1s=1;           % スレーブ群の最良位置に対する重み 
c2s=1;           % 全体の最良位置に対する重み
c1m=1;           % マスター群の最良位置に対する重み
c2m=1;           % スレーブ群の最良位置に対する重み
c3m=1;           % 全体の最良位置に対する重み

maxT=500;     % 繰り返し回数
D=2;             % 次元

%****************　[STEP1] 初期粒子群の発生 ********************************
 
%スレーブ粒子
sp = 0.5-rand(ns,2);  %初期粒子群発生 
sp(:,3) = 0;
sp_v = rand(ns,2);  
sp_best = sp; %ローカルベスト   

%マスター粒子
mp = 0.5-rand(nm,2);  %初期粒子群発生 
mp(:,3) = 0;
mp_v = rand(nm,2);  
mp_best = mp; %ローカルベスト   

P_S_G(1,D+1) = 0; %スレーブ粒子群全体のグローバルベスト
P_G(1,D+1) = 0; %粒子群全体のグローバルベスト

f0 = figure;
plot(sp(:,1),sp(:,2),'o')
hold on
plot(mp(:,1),mp(:,2),'o')
hold off
xlim([-0.5 0.5])
ylim([-0.5 0.5])

for count=1:maxT
%****************　[STEP2] 評価値計算 **************************************
            
    for t=1:ns
        sp(t,3) = 1/abs(sp(t,1))+abs(sp(t,2));
    end
    for t=1:ns
        mp(t,3) = 1/abs(mp(t,1))+abs(mp(t,2));
    end    
               
%****************　[STEP3] 各粒子の最良位置の更新 ***************************
            
    for i=1:ns 
        if(sp(i,3) > sp_best(i,3))
            for j=1:D+1
               sp_best(i,j) = sp(i,j);
            end
        end
    end 
    
    for i=1:ns 
        if(mp(i,3) > mp_best(i,3))
            for j=1:D+1
               mp_best(i,j) = mp(i,j);
            end
        end
    end   
           
 %****************　[STEP4] 全体の最良位置の更新 ***************************
    for i=1:ns %スレーブ群全体
        if(sp(i,3) > P_S_G(1,3))
            for j=1:D+1
                P_S_G(1,j) = sp(i,j);
            end
        end 
    end
    
    for i=1:nm %全体の最良位置
        if(mp(i,3) > P_G(1,3))
            for j=1:D+1
                P_G(1,j) = mp(i,j);
            end
         end 
    end
    
    if(P_S_G(1,3) > P_G(1,3))
        for j=1:D+1
            P_G(1,j) = P_S_G(1,j);
        end
    end        
                 
 %****************　[STEP5] 各粒子の速度，位置を計算 ***************************     
            
    %パラメータ更新
    for i=1:nm
        for j=1:D
            r1 = rand/100;
            r2 = rand/100;
            r3 = rand/100;
            mp_v(i,j) = w*mp_v(i,j) + c1m*r1*(mp_best(i,j)-mp(i,j)) + c2m*r2*(P_S_G(1,j)-mp(i,j)) + c3m*r3*(P_G(1,j) - mp(i,j)); %速度更新
            mp(i,j) = mp(i,j) + mp_v(i,j); %位置更新
        end
    end
    for i=1:ns
        for j=1:D
            r1 = rand/100;
            r2 = rand/100;
            sp_v(i,j) = c1s*r1*(sp_best(i,j)-sp(i,j)) + c2s*r2*(P_G(1,j)-sp(i,j)); %速度更新
            sp(i,j) = sp(i,j) + sp_v(i,j); %位置更新
        end
    end
    

    if (mod(count,250)==0) %250更新ごとに図を表示
        f1 = figure;
        plot(sp(:,1),sp(:,2),'o')
        hold on
        plot(mp(:,1),mp(:,2),'o')
        hold off
        xlim([-0.5 0.5])
        ylim([-0.5 0.5])
    end
       
end

Tags:

群知能