一、粒子群算法簡介
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算法的流程和僞代碼
二、源代碼
clear
clc
close all
%% 參數初始化
c1=2.05;
c2=2.05;
maxgen=5000;
sizepop=30;
k=0.6;
% wV=1.1;
wP=1.1;
v=5;
popmax=30;
popmin=-30;
% pso_option = struct('c1',1.5,'c2',1.7,'maxgen',200,'sizepop',20, ...
% 'k',0.6,'wV',1,'wP',1,'v',5, ...
% 'popcmax',100,'popcmin',0.1,'popgmax',10^3,'popgmin',10^(-2));
% c1:初始為1.5,pso參數局部搜尋能力
% c2:初始為1.7,pso參數全局搜尋能力
% maxgen:初始為200,最大進化數量
% sizepop:初始為20,種群最大數量
% k:初始為0.6(k belongs to [0.1,1.0]),速率和x的關系(V = kX)
% wV:初始為1(wV best belongs to [0.8,1.2]),速率更新公式中速度前面的彈性系數
% wP:初始為1,種群更新公式中速度前面的彈性系數
% v:初始為5,SVM Cross Validation參數
% popcmax:初始為100,SVM 參數c的變化的最大值.
% popcmin:初始為0.1,SVM 參數c的變化的最小值.
% popgmax:初始為1000,SVM 參數g的變化的最大值.
% popgmin:初始為0.01,SVM 參數c的變化的最小值.
D=10; %%%維數
Vmax =k*popmax;
Vmin = -Vmax ;
eps =1E-5;
%% 産生初始粒子和速度
pop=zeros(sizepop,D);
V=zeros(sizepop,D);
fitness=zeros(sizepop,1);
for i=1:sizepop
% 随機産生種群和速度
pop(i,:) = (popmax-popmin)*rand(1,D)+popmin;
V(i,:)=Vmax*rands(1,D);
% 計算初始适應度
fitness(i)=myfunc_fit1(pop(i,:));
end
Xd_ave0=repmat(sum(pop)/sizepop,sizepop,1);
D_t0=sum((sum((pop-Xd_ave0).^2,2)).^0.5)/sizepop/(popmax-popmin);
wV=1/(1+exp(-12*(D_t0-0.5)));
D_min=D_t0*0.2;
% 找極值和極值點
[global_fitness bestindex]=min(fitness); % 全局極值
local_fitness=fitness; % 個體極值初始化
global_x=pop(bestindex,:); % 全局極值點
local_x=pop; % 個體極值點初始化
% 每一代種群的平均适應度
avgfitness_gen = zeros(maxgen,1);
fit_gen=zeros(maxgen,1);
%% 疊代尋優
for i=1:maxgen
for j=1:sizepop
%速度更新
V(j,:) = wV*V(j,:) + c1*rand*(local_x(j,:) - pop(j,:)) + c2*rand*(global_x - pop(j,:));
if find(V(j,:) > Vmax)
V_maxflag=find(V(j,:) > Vmax);
V(j,V_maxflag) = Vmax;
end
if find(V(j,1) < Vmin)
V_minflag=find(V(j,1) < Vmin);
V(j,V_minflag) = Vmin;
end
%種群更新
pop(j,:)=pop(j,:) + wP*V(j,:);
if find(pop(j,:) > popmax)
pop_maxflag=find(pop(j,:) > popmax);
pop(j,pop_maxflag) = popmax;
end
if find(pop(j,:) < popmin)
pop_minflag=find(pop(j,:) < popmin);
pop(j,pop_minflag) = popmin;
end
% 自适應粒子變異
if rand>0.5
k=ceil(2*rand);
pop(j,k) = (popmax-popmin)*rand + popmin;
end
%适應度值
fitness(j)=myfunc_fit1(pop(j,:));
%個體最優更新
if fitness(j) < local_fitness(j)
local_x(j,:) = pop(j,:);
local_fitness(j) = fitness(j);
end
% if fitness(j) == local_fitness(j) && length(pop(j,:) < local_x(j,:))
% local_flag=find(pop(j,:) < local_x(j,:));
% local_x(j,local_flag) = pop(j,local_flag);
% local_fitness(j) = fitness(j);
% end
%群體最優更新
if fitness(j) < global_fitness
global_x = pop(j,:);
global_fitness = fitness(j);
end
function C=PSO_FUNC(X)
global G_AC T_C hour_num Wind_V
C_W=110; %%%風力發電
P_W=100;
u_PW=6;
m_WG=20;
r0=0.06;
v_ci=3; %切入風速
v_r=12; %額定風速
v_co=25; %切除風速
P_r=P_W;
V_t=repmat(Wind_V(:,3),52,1);
C_S=0.7; %%%%光伏發電
P_S=0.2;
u_PS=0.009;
m_PV=25;
C_B=0.5; %%%蓄電池
u_WB=0.0014;
m_B=5;
sigam_bat=1e-4; %%%自放電率
N_B=2000;
W_B=0.64;
Wbat_0=0.5*N_B*W_B;
Pbat_max=0.2*N_B*W_B;
Pbat_min=-0.2*N_B*W_B;
Pbat_maxt=Pbat_max;
Pbat_mint=Pbat_min;
Wbat_t=zeros(hour_num,1);
Wbat_t(1)=Wbat_0; %%%剩餘電量
Pbat_t=zeros(hour_num,1);
Pbat_t(1)=Pbat_max;
C_d=10; %%%柴油機發電
u_Pd=0.95;
P=4.62;
Q_d0=0.22;
m_die=10;
P_STC=0.2;
G_STC=1;
K=-0.47;
Tr=298.15;
T_C=T_C+273.15;
P1_t=300; %%%居民負荷峰值 KW
Pdes_t=200; %%%海水淡化負荷 KW
P_des=25; %%%單台海水淡化機組的額定功率 KW
N_des=8; %%%海水淡化機組總台數
G_des=100/24; %%%單台機組每小時的産水量 100t/d
Rwater_t=500/24; %%%島上全天用水需求 t
Rdes_min=0;
Rdes_max=8*100/24;
eta_c=0.97;
Rdes_t=zeros(hour_num,1);
Rdes_t(1)=Rdes_max; %%最初蓄水量
P_PV=zeros(hour_num,1);
P_WG=zeros(hour_num,1);
P_PVM_t=zeros(hour_num,1);
P_WGM_t=zeros(hour_num,1);
P_net_t=zeros(hour_num,1);
Pdie_t=zeros(hour_num,1);
C_f=0; %%%柴油年成本
yeushu1=0;
yeushu2=0;
yeushu3=0;
yeushu4=0;
yeushu5=0;
yeushu6=0;
yeushu7=0;
yeushu8=0;
for i=1:hour_num
P_PV(i)=P_STC*G_AC(i).*(1+K*(T_C(i)-Tr))/G_STC;
a=P_r/(v_r^3-v_ci^3);
b=v_ci^3/(v_r^3-v_ci^3);
if (V_t(i)<v_ci)
P_WG(i)=0;
elseif (v_ci<V_t(i)<v_r)
P_WG(i)=a*V_t(i)^3-b*P_r;
elseif (v_r<V_t(i)<v_co)
P_WG(i)=P_r;
else
P_WG(i)=0;
end
if(Rdes_t(i)-Rdes_max>=Rwater_t)
Ndes_mint=0;
else
Ndes_mint=(Rwater_t-(Rdes_t(i)-Rdes_max))/G_des;
end
if(Rdes_t(i)+N_des*P_des-Rwater_t<=Rdes_max)
Ndes_maxt=N_des*P_des;
else
Ndes_maxt=(Rdes_max+Rwater_t-Rdes_t(i))/P_des;
end
Pdes_mint=Ndes_mint*P_des;
Pdes_maxt=Ndes_maxt*P_des;
P_PVM_t(i)=P_PV(i);
P_WGM_t(i)=P_WG(i);
P_net_t(i)=P_WGM_t(i)+P_PVM_t(i)-P1_t;
if(i>=2)
Pbat_maxt=min([Pbat_max ((1-sigam_bat)*Wbat_t(i)-((1-sigam_bat)*Wbat_t(i)+Pbat_t(i)))*eta_c]);
Pbat_mint=max([Pbat_min ((1-sigam_bat)*Wbat_t(i)-((1-sigam_bat)*Wbat_t(i)+Pbat_t(i)))/eta_c]);
end
三、運作結果
四、matlab版本及參考文獻
1 matlab版本
2014a
2 參考文獻
[1] 包子陽,餘繼周,楊杉.智能優化算法及其MATLAB執行個體（第2版）[M].電子工業出版社，2016.
[2]張岩,吳水根.MATLAB優化算法源代碼[M].清華大學出版社，2017.