一、簡介
粒子群演算法源於複雜適應系統（Complex Adaptive System,CAS）。CAS理論於1994年正式提出，CAS中的成員稱為主體。比如研究鳥群系統，每個鳥在這個系統中就稱為主體。主體有適應性，它能夠與環境及其他的主體進行交流，並且根據交流的過程“學習”或“積累經驗”改變自身結構與行為。整個系統的演變或進化包括：新層次的產生（小鳥的出生）；分化和多樣性的出現（鳥群中的鳥分成許多小的群）；新的主題的出現（鳥尋找食物過程中，不斷髮現新的食物）。
所以CAS系統中的主體具有4個基本特點（這些特點是粒子群演算法發展變化的依據）：
首先，主體是主動的、活動的。
主體與環境及其他主體是相互影響、相互作用的，這種影響是系統發展變化的主要動力。
環境的影響是巨集觀的，主體之間的影響是微觀的，巨集觀與微觀要有機結合。
最後，整個系統可能還要受一些隨機因素的影響。
粒子群演算法就是對一個CAS系統－－－鳥群社會系統的研究得出的。
粒子群演算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy於1995年提出，它的基本概念源於對鳥群覓食行為的研究。設想這樣一個場景:一群鳥在隨機搜尋食物，在這個區域裡只有一塊食物，所有的鳥都不知道食物在哪裡，但是它們知道當前的位置離食物還有多遠。那麼找到食物的最優策略是什麼呢?最簡單有效的就是搜尋目前離食物最近的鳥的周圍區域。
PSO演算法就從這種生物種群行為特性中得到啟發並用於求解優化問題。在PSO中，每個優化問題的潛在解都可以想象成d維搜尋空間上的一個點，我們稱之為“粒子”（Particle），所有的粒子都有一個被目標函式決定的適應值(Fitness Value )，每個粒子還有一個速度決定他們飛翔的方向和距離，然後粒子們就追隨當前的最優粒子在解空間中搜索。Reynolds對鳥群飛行的研究發現。鳥僅僅是追蹤它有限數量的鄰居但最終的整體結果是整個鳥群好像在一箇中心的控制之下.即複雜的全域性行為是由簡單規則的相互作用引起的。
二、原始碼

clc,clear,close all
warning off
format longG
tic;
global Nseat Dab Dbc Dac fab fbc fac ti ta tb
% 載入引數
parameters;    % 系統的引數
% 第index_t期
index_t = 10; 
Vt_last = 0;
% PSO 引數
c1 = 1.4995;  
c2 = 1.4995;
Vmin = -1;
Vmax = 1;
maxiter = 30;  % 迭代次數
sizepop = 50;  % 種群數量
nvar = 3;      % 3個變數
popmin = 0;  popmax = max([ta(index_t),tb(index_t)]); % x
% 初始化種群
for i=1:sizepop
    pop(i,:) = popmin + (popmax-popmin).*rand(1,nvar);
    fitness(i) = fun(pop(i,:), Vt_last, index_t);
    V(i,1:3) = 0;
end
% 記錄一組最優值
[bestfitness,bestindex]=min(fitness);
zbest=pop(bestindex,:);   % 全域性最佳
gbest=pop;                % 個體最佳
fitnessgbest=fitness;     % 個體最佳適應度值
fitnesszbest=bestfitness; % 全域性最佳適應度值
wmax = 0.9;  wmin = 0.4;
% 迭代尋優
for i=1:maxiter
    for j=1:sizepop
        % 自適應權重1
%         w = wmin + (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( max(fitness)-min(fitness) );
        % 自適應權重2
%         w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
        % 自適應權重3
        if fitnessgbest(j)<=mean(fitness)
            w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
        else
            w = wmax;
        end
        % 速度更新
        V(j,:) = w*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        % V--防止越界
        for k=1:nvar
            if V(j,k)>Vmax
                V(j,k)=Vmax;
            end
            if V(j,k)<Vmin
                V(j,k)=Vmin;
            end
        end

        % 個體更新
        pop(j,:) = pop(j,:) + 0.5 * V(j,:);
        
        % x  越界限制
        for k=1:nvar
            if pop(j,k)>popmax
                pop(j,k)=popmax;
            end
            if pop(j,k)<popmin
                pop(j,k)=popmin;
            end
        end
        
        % 適應度更新
        fitness(j) = fun(pop(j,:), Vt_last, index_t);
        
        % 比較  個體間比較
        if fitness(j)<fitnessgbest(j)
            fitnessgbest(j) = fitness(j);
            gbest(j,:) = pop(j,:);
        end
        if fitness(j)<bestfitness
            bestfitness = fitness(j);
            zbest =  pop(j,:);
        end
        
    end
    fitness_iter(i) = bestfitness;
    
end
function [zbest,fitness_iter,bestfitness] = APSO_Vt(index_t,Vt_last, maxiter,sizepop  )
global Nseat Dab Dbc Dac fab fbc fac ti ta tb
% 載入引數
% parameters;    % 系統的引數
% 第index_t期
% index_t = 10; 
% Vt_last = 0;
% PSO 引數
c1 = 1.4995;  
c2 = 1.4995;
Vmin = -1;
Vmax = 1;
% maxiter = 20;  % 迭代次數
% sizepop = 20;  % 種群數量
nvar = 3;      % 3個變數
popmin = 0;  popmax = max([ta(index_t),tb(index_t)]); % x
% 初始化種群
for i=1:sizepop
    pop(i,:) = fix(popmin + (popmax-popmin).*rand(1,nvar));
    fitness(i) = fun(pop(i,:), Vt_last, index_t);
    V(i,1:3) = 0;
end
% 記錄一組最優值
[bestfitness,bestindex]=min(fitness);
zbest=pop(bestindex,:);   % 全域性最佳
gbest=pop;                % 個體最佳
fitnessgbest=fitness;     % 個體最佳適應度值
fitnesszbest=bestfitness; % 全域性最佳適應度值
wmax = 0.9;  wmin = 0.4;
% 迭代尋優
for i=1:maxiter
    for j=1:sizepop
        % 自適應權重1
%         w = wmin + (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( max(fitness)-min(fitness) );
        % 自適應權重2
%         w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
        % 自適應權重3
        if fitnessgbest(j)<=mean(fitness)
            w = wmin - (wmax-wmin)*(fitnessgbest(j)-min(fitness))/( mean(fitness)-min(fitness) );
        else
            w = wmax;
        end
        % 速度更新
        V(j,:) = w*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        % V--防止越界
        for k=1:nvar
            if V(j,k)>Vmax
                V(j,k)=Vmax;
            end
            if V(j,k)<Vmin
                V(j,k)=Vmin;
            end
        end

        % 個體更新
        pop(j,:) = fix(pop(j,:) + 0.5 * V(j,:));
        
        % x  越界限制
        for k=1:nvar
            if pop(j,k)>popmax
                pop(j,k)=popmax;
            end
            if pop(j,k)<popmin
                pop(j,k)=popmin;
            end
        end
        function y = fun(x, Vt_last, index_t )
global Nseat Dab Dbc Dac fab fbc fac ti ta tb

xab = x(1);
xbc = x(2);
xac = x(3);

gt = fab*min(Dab(index_t),xab) + fbc*min(Dbc(index_t),xbc) + fac*min(Dac(index_t),xac);
Vt = gt + Vt_last;   % 目標值越大越好
if( (xab+xac<=ta(index_t)) && (xbc+xac<=tb(index_t)) )
    y = 1./Vt;           % 求最小值，y越小越好
else
    y = 1e6;             % 不滿足約束條件，則y為極大值
end


 

三、執行結果

四、備註
版本：2014a