粒子群最佳化演算法在複雜函數參數估計中的應用

全文：http://www.lingch.net/db/download.asp?tab=softdown&fild=9&id=30 摘要：粒子群最佳化（Particle Swarm Optimization, PSO）演算法是有效現代啟發學習法搜尋演算法之一,其本身有計算簡便、收斂速度快和收斂準確等特點。本文將其應用於複雜函數的參數估計中，利用本文提及的基於PSO演算法複雜函數參數估計法，可以準確估計複雜函數的參數，且具有計算簡便、收斂速度快、估計準確和高適用性等特點。 關鍵字：粒子群 ; 最佳化 ; 啟發學習法演算法 ; PSO ; 參數估計 附錄 MatLab來源程式
%------------------------------主程式
%粒子群最佳化演算法在參數最小二乘估計中的應用
%陳毓靈
%Post2Ling@hotmail.com

points=50; %樣本點數
vec=4; %待估參數向量維數

Pop=200; %種群大小
Vmax=20; %最大速度
Vmin=-20; %最小速度
nIter=150; %迭代次數
Wmax=0.8; %最大慣性因子
Wmin=0.0; %最小慣性因子
Wn=1.0; %慣性因子衰減指數
FaiI=0.3; %自身學習率
FaiG=0.8; %社會學習率

%輸入樣本 X Y
X=51+rand(50,1)*50;
Y = 100+(200./(1+exp(0.3.*(X-75))))+(rand(50,1)-0.5)*20;

%粒子群
%當前位置vec維 個體曆史最優位置vec維 全域最優位置vec維 當前值 個體曆史最優值 全域最優值
P=rand(Pop,vec*3+3)*300;
P(:,vec+1:vec*2)=P(:,1:vec);

% 更新適應值
P(:,vec*3+1)=Fitness(P(:,1:vec),X,Y);
P(:,vec*3+2)=P(:,vec*3+1);
P(:,vec*3+3)=P(1,vec*3+1);
% 更新全域最優位置和適應值
for j=1:Pop
if P(j,vec*3+2)<P(1,vec*3+3)
P(1,vec*3+3)=P(j,vec*3+2);
P(1,vec*2+1:vec*3)=P(j,vec+1:vec*2);
end
end
P(:,vec*2+1:vec*3)=ones(Pop,1)*P(1,vec*2+1:vec*3);
P(:,vec*3+3)=P(1,vec*3+3);

%粒子群速度
V=(rand(Pop,vec)-0.5)*50;

for i=1:nIter

% 更新慣性因子
W=((nIter-i)^Wn)/(nIter^Wn)*(Wmax-Wmin)+Wmin;

% 更新速度
V=W.*V + FaiI.*(P(:,vec+1:vec*2)-P(:,1:vec))+FaiG.*(P(:,vec*2+1:vec*3)-P(:,1:vec));
% 限制速度
V(find(V>Vmax))=Vmax;
V(find(V<Vmin))=Vmin;

% 更新位置
P(:,1:vec)=P(:,1:vec) + V(:,1:vec);

% 更新適應值
P(:,vec*3+1)=Fitness(P(:,1:vec),X,Y);

% 更新個體最優位置和適應值
for j=1:Pop
if P(j,vec*3+1) < P(j,vec*3+2)
P(j,vec+1:vec*2)=P(j,1:vec);
P(j,vec*3+2)=P(j,vec*3+1);
end
end

% 更新全域最優位置和適應值
for j=1:Pop
if P(j,vec*3+2)<P(1,vec*3+3)
P(1,vec*3+3)=P(j,vec*3+2);
P(1,vec*2+1:vec*3)=P(j,vec+1:vec*2);
end
end
% 複製冗餘部分以方便計算
P(:,vec*2+1:vec*3)=ones(Pop,1)*P(1,vec*2+1:vec*3);
P(:,vec*3+3)=P(1,vec*3+3);

end

%輸出結果
% [P(1,vec*2+1:vec*3),P(1,vec*3+3)]
[P(1,vec*2+1),P(1,vec*2+2),P(1,vec*2+3)/300,P(1,vec*2+4)]

BP=51:100;
HR = P(1,vec*2+1)+(P(1,vec*2+2)./(1+exp((P(1,vec*2+3)./300).*(BP-P(1,vec*2+4)))));%+(rand(1,50)-0.5)*10;

h=plot(X,Y,'r*',BP,HR,'b');
xlabel('x')
ylabel('Logistic')
title('4參數Logistic函數映像與樣本比較 (p1=100;p2=200;p3=0.3;p4=75)')
legend('樣本','估計函數')

%----------------------適應值函數
function y = Fitness( P,X,Y)
%最小二乘估計適應值函數
%以誤差二乘方和作為適應值，目標是最小化適應值

p_size=size(P);
x_size=size(X);
y_size=size(Y);

%結果
y=zeros(p_size(1),1);

for i=1:p_size(1)%對每個粒子
for j=1:x_size(1)%對每個樣本點
y(i)=y(i)+((Y(j)-(P(i,1)+(P(i,2)./(1+exp((P(i,3)./300).*(X(j,1)-P(i,4)))))))^2);%y(i)為i例子得當前適應值（誤差二乘方和）
end
end

% X(:,p_size(2))=1;%擴充矩陣
%sum(P(i,:).*X(j,:))