1.演算法概述

        粒子群優化演算法(PSO)是一種進化計算技術(evolutionary computation)，1995 年由Eberhart 博士和kennedy 博士提出，源於對鳥群捕食的行為研究 。該演算法最初是受到飛鳥叢集活動的規律性啟發，進而利用群體智慧建立的一個簡化模型。粒子群演算法在對動物叢集活動行為觀察基礎上，利用群體中的個體對資訊的共享使整個群體的運動在問題求解空間中產生從無序到有序的演化過程，從而獲得最優解。  

       非線性系統的辨識一直是現代辨識領域中的一個主要課題。針對非線性系統中Hammerstein 模型, 討論了PSO 的基本演算法與引數初值的設定與選擇方法。通過模擬實驗說明

PSO 過程:

第1 步: 種群隨機初始化。

第2 步: 對種群內的每一個個體計算適應值(fitness value) , 適應值與最優解的距離直接有關。

第3 步: 種群根據適應值進行復制。

第4 步: 如果終止條件滿足, 則停止; 否則轉到第2 步。

       PSO 演算法中每個優化問題的解都有是搜尋空間中的一隻鳥, 稱為粒子。與其他進化計算技術不同的是群體中的每個粒子可以記憶自己到過的最優位置, 並能感知鄰近群體已達到的最優位置, 每個粒子能夠根據自身到過的最優位置和鄰近群體已到過的最優位置來更新自己

       粒子的速度和位置根據以下的公式從h 代更新到h+ 1代:

 Hammerstein模型

       將Hammerstein模型轉換為一類中間模型。然後,提出利用一種改進的粒子群優化（improved particle swarm optimization,IPSO）演算法獲得中間模型的引數估計值。接著,通過相應的數學關係來達到對Hammerstein模型的辨識。

2.模擬效果預覽

3.核心MATLAB程式

function [pls_P,pls_Q,pls_W,pls_B]=pls_model(X_M, Y_M, iter)
 
% X=TP'+E=∑tp'+E
% Y=UQ'+F=∑uq'+F
% u=bt[pls_P,pls_Q,pls_W,pls_B]=pls_model(X_M, Y_M, iter)
% b=u't/t't
% Y=TBQ'+F
 
% pls_P the loadings for every pls component of X
% pls_Q the loadings for every pls componentof Y
% pls_W the weights
% pls_B the regression coefficients for the inner relation.
[N,K]=size(X_M);
[N,P]=size(Y_M);
 
if iter > min(N,K)
    iter = min(N,K);
end
 
% start algrithm: PLS
for counter = 1:iter; % extract 3 PLS components
    pls_u = Y_M(:,1); % get any column from Y
    t_old_n = 1.0;% the ini of t old
    t_new_n = 0.0;% the ini of t new
    pls_t = 0;
    while ( abs( t_old_n - t_new_n )>0.0001)
        pls_t_old = pls_t;
        % in the X block
        pls_w = (pls_u' * X_M / ( pls_u' * pls_u ))';
        pls_w = pls_w/norm(pls_w);
        pls_t = X_M * pls_w;
        % in the Y block
        if(P>1)
            pls_q = (pls_t' * Y_M /(pls_t'*pls_t))';
            pls_q = pls_q/norm(pls_q);
            pls_u = Y_M*pls_q;
        else
            pls_q=1;
        end
        % compare t,check for convergence
        t_old_n = norm( pls_t_old );
        t_new_n = norm( pls_t);
    end
 
    pls_p = (pls_t'*X_M/(pls_t'*pls_t))';
    pls_p_n=norm(pls_p);    
    pls_p = pls_p/pls_p_n;
    pls_t = pls_t*pls_p_n;
    pls_w = pls_w*pls_p_n;
    pls_b = pls_u'*pls_t/(pls_t'*pls_t);
 
    X_M = X_M - pls_t*pls_p';
    Y_M = Y_M - pls_b*pls_t*pls_q';
    pls_P(:,counter)=pls_p;
    pls_Q(:,counter)=pls_q;
    pls_W(:,counter)=pls_w;
    pls_B(:,counter)=pls_b;
end
02-001m