一、粒子群算法簡介
1 粒子群算法的概念
粒子群優化算法(PSO：Particle swarm optimization) 是一種進化計算技術（evolutionary computation）。源于對鳥群捕食的行為研究。粒子群優化算法的基本思想：是通過群體中個體之間的協作和資訊共享來尋找最優解．
PSO的優勢：在于簡單容易實作并且沒有許多參數的調節。目前已被廣泛應用于函數優化、神經網絡訓練、模糊系統控制以及其他遺傳算法的應用領域。
2 粒子群算法分析
2.1基本思想
粒子群算法通過設計一種無品質的粒子來模拟鳥群中的鳥，粒子僅具有兩個屬性：速度和位置，速度代表移動的快慢，位置代表移動的方向。每個粒子在搜尋空間中單獨的搜尋最優解，并将其記為目前個體極值，并将個體極值與整個粒子群裡的其他粒子共享，找到最優的那個個體極值作為整個粒子群的目前全局最優解，粒子群中的所有粒子根據自己找到的目前個體極值和整個粒子群共享的目前全局最優解來調整自己的速度和位置。下面的動圖很形象地展示了PSO算法的過程：
2 更新規則
PSO初始化為一群随機粒子(随機解)。然後通過疊代找到最優解。在每一次的疊代中，粒子通過跟蹤兩個“極值”(pbest，gbest)來更新自己。在找到這兩個最優值後，粒子通過下面的公式來更新自己的速度和位置。
公式(1)的第一部分稱為【記憶項】，表示上次速度大小和方向的影響；公式(1)的第二部分稱為【自身認知項】，是從目前點指向粒子自身最好點的一個矢量，表示粒子的動作來源于自己經驗的部分；公式(1)的第三部分稱為【群體認知項】，是一個從目前點指向種群最好點的矢量，反映了粒子間的協同合作和知識共享。粒子就是通過自己的經驗和同伴中最好的經驗來決定下一步的運動。以上面兩個公式為基礎，形成了PSO的标準形式。
公式(2)和 公式(3)被視為标準PSO算法。
3 PSO算法的流程和僞代碼
二、部分源代碼
%------主函數源程式－周排程（main.m）
%------基本粒子群優化算法（Particle Swarm Optimization）-----------
%------名稱：基本粒子群優化算法（PSO）
%------作用：求解優化問題(周排程)
%------說明：全局性，并行性，高效的群體智能算法
%------初始格式化--------------------------------------------------
tic
clear all;
clc;
format short;
%% ------給定水電站初始化條件----------------------------------------------
%------------水電站1---------------------------------------
Vmax1=9015*10^4; %水庫容量上限(m3)
Vmin1=7000*10^4; %水庫容量下限(m3)
H1=640; %水庫容量初始值水位(m)
V1=(2.0554*H1^2-2413.5002*H1+709934.65)*10^4 ; %水庫庫容與水位的關系
h1=91; %初始水庫水頭(m)
qr(:,1:7)=[20.4 25.2 22.1 19.3 16.4 23.3 28.6 ]; %水庫來水流量(m3/s)
qmax1=44; %水庫引用流量上限(m3/s)
qmin1=0; %水庫引用流量下限(m3/s)
A1=9.8*10^3; %水庫出力系數
k1=0.65; %發電效率
t=8.64*10^4; %水庫發電引用流量時間段(s)
%-----------水電站2--------------------------------------
Vmax2=2080*10^4; %水庫容量上限(m3)
Vmin2=1530*10^4; %水庫容量下限(m3)
H2=540; %水庫容量初始值水位(m)
V2=(2.6408*H2^2-2763.6392*H2+724014.90)*10^4; %水庫庫容與水位的關系
h2=57; %初始水庫水頭(m)
qr(:,8:14)=[22.4 18.3 26.4 25.2 17.6 24.6 27.2]; %水庫來水流量(m3/s)
qmax2=32; %水庫引用流量上限(m3/s)
qmin2=0; %水庫引用流量下限(m3/s)
A2=9.8*10^3; %水庫出力系數
k2=0.6; %發電效率
%% 給定粒子群計算初始化條件
T=7; %周期
c1=2; %學習因子1
c2=2; %學習因子2
w1=0.9; %慣性權重
w2=0.4;
MaxIter=100; %最大疊代次數
Iter=0;
vmax=2; %最大運動速度
vmin=0; %最小運動速度
N=20; %初始化群體個體數目
%------初始化種群的個體(可以在這裡限定位置和速度的範圍)------------
q=zeros(N,2*T); %位置随機初始化位置
v=rand(N,2*T); %速度随機初始化速度
%------适應度函數源程式（fitness.m）----------------------------
Fx1=sum((k1*A1*q(:,1:7).*h1*t/(3.6*10^10))')'; %水庫1總發電量（萬kWh）
Fx2=sum((k2*A2*q(:,8:14).*h2*t/(3.6*10^10))')'; %水庫2總發電量（萬kWh）
Fx=Fx1+Fx2; %水庫總發電量（萬kWh）
%------設定目前粒子的最好位置----------------------------------
Pbest=q;
FPbest=Fx;
%------找到初始粒子的最好粒子------------------------------------
[Fgbest,r]=max(Fx);
E=Fgbest;%記錄目前全局最優值
Best=q(r,:);%用于儲存最優粒子位置
%% ------------------------------循環-----------------------------
while(Iter<=MaxIter)
Iter=Iter+1;
%------------------------------更新慣性權重和速度----------------------
w=((w1-w2)*(MaxIter-Iter)/MaxIter)+w2;% 典型的線性遞減政策
% 線性微分遞減政策
% 先增後減政策
% 非線性慣性權重政策
R1=rand(N,2*T);
R2=rand(N,2*T);
%-------------------------粒子速度更新----------------------------
B=repmat(q(r,:),N,1);
v=w*v+c1*R1.*(Pbest-q)+c2*R2.*(B-q);
%-------------------------粒子速度限制----------------------------
changev=v>vmax;
v(find(changev))=vmax;
changev=v<vmin;
v(find(changev))=vmin;
%-------------------------粒子流量更新----------------------------
q=q+1.0*v;
%-------------------------各時段引用流量限制----------------------------
for i=1:7
q1=q(:,i);
changeq1=q1>qmax1;
q1(find(changeq1))=qmax1;
changeq1=q1<qmin1;
q1(find(changeq1))=qmin1;
q(:,i)=q1;
end
for i=8:14
q2=q(:,i);
changeq2=q2>qmax1;
q2(find(changeq2))=qmax1;
changeq2=q2<qmin1;
q2(find(changeq2))=qmin1;
q(:,i)=q2;
end
%-------------------------各時段庫容大小(m3)----------------------------
qr1=qr(:,1:7);
qr2=qr(:,8:14);
Qr1=repmat(qr1,N,1); %----各時段來水大小
Vs1=repmat(0,N,T+1); %----各時段初始庫容大小
Vj1=repmat(0,N,T); %----各時段庫容大小(首末平均值)
quit1=repmat(0,N,T); %----各時段棄水大小
Qr2=repmat(qr2,N,1); %----各時段來水大小
Vs2=repmat(0,N,T+1); %----各時段初始庫容大小
Vj2=repmat(0,N,T); %----各時段庫容大小(首末平均值)
quit2=repmat(0,N,T); %----各時段棄水大小
%-------------------------各時段水頭大小----------------------------------
for i=1:N
Vs1(i,1)=V1;
Vs2(i,1)=V2;
end
for i=1:N
for j=2:(T+1)
q1=q(:,1:7);
Vs1(i,j)=Vs1(i,(j-1))+(Qr1(i,(j-1))-q1(i,(j-1)))*t;
if Vs1(i,j)>=Vmax1
quit1(i,j-1)= Vs1(i,j)-Vmax1;
Vs1(i,j)=Vmax1;
end
if Vs1(i,j)<=Vmin1
q1(i,j-1)=(Vs1(i,j)-Vmin1+Qr1(i,j-1)*t)/(t);
Vs1(i,j)=Vmin1;
end
end
end
for i=1:N
for j=2:(T+1)
q2=q(:,8:14);
Vs2(i,j)=Vs2(i,(j-1))+(Qr2(i,(j-1))-q2(i,(j-1)))*t+quit1(i,(j-1))+q1(i,j-1)*t;
if Vs2(i,j)>=Vmax2
quit2(i,j-1)= Vs2(i,j)-Vmax2;
Vs2(i,j)=Vmax2;
end
if Vs2(i,j)<=Vmin2
q2(i,j-1)=(Vs2(i,j)-Vmin2+Qr2(i,j-1)*t)/(t);
Vs2(i,j)=Vmin2;
end
end
end
for i=1:N
for j=1:T
Vj1(i,j)=(Vs1(i,j)+Vs1(i,j+1))/2;
Vj2(i,j)=(Vs2(i,j)+Vs2(i,j+1))/2;
end
end