标準PSO辨識NARMAX模型源碼程式

% 用标準PSO辨識NARMAX模型

% y=[-0.4 0.2 0.4 0.8 0.2 0.3]*...

%   [y0(k-1) y0(k-2)*y0(k-3) y0(k-4) u(k-1)^3 y0(k-2)^2 u(k-3)]'

% +e(k)

% function [iter,Xgbest,fgbest]=sPso(err,var,Nc_max)

clc,clear,format long

%-----------------------------------------------------------paramters setup

var=0.001;         % 噪聲方差

% iter_max=Nc_max*4;

err=0.01;

iter_max=1000;    % 最大疊代次數

L=15;           % 視窗寬度

N=30;            % 種群規模

C1=2;            % 加速度常數

C2=C1;

Xmin=-3;         % 解取值範圍[Xmin,Xmax]

Xmax=3;

p=6;             % 粒子維數

w=linspace(0.9,0.5,iter_max);   % 慣性權重

X=Xmin+(Xmax-Xmin)*rand(p,N);   % 粒子位置

Xpbest=X;                       % 個體最佳位置

Xgbest=Xmin+(Xmax-Xmin)*rand(p,1);   % 種群最佳位置

fpbest=0*rand(1,N);           % 個體最佳适應度值

fgbest=0;                     % 種群最佳适應度值

fgbest_fig=zeros(1,iter_max); 

Xgbest_fig=zeros(p,iter_max);

Vmax=(Xmax-Xmin)*0.2;

V=Vmax*(2*rand(p,N)-1);

u=idinput(L,'rgs',[0 1],[-1 1]);      % 随機白噪聲序列，取L個，均值為0，方差為1

e=idinput(L,'rgs',[0 1],[-var var]);  % 高斯白噪聲，均值為0，方差為var

theta0=[-0.4 0.2 0.4 0.8 0.2 0.3];    % 待辨識參數真值

%-----------------------------------------------output of theoretical value 

y0(1:4)=0;

y(1:4)=0;

for k=5:L

    y0(k)=theta0*[y0(k-1) y0(k-2)*y0(k-3) y0(k-4) u(k-1)^3 y0(k-2)^2 u(k-3)]'+e(k);

end

%----------------------------------------------------------------------main

iter=0;

while iter<iter_max

    iter=iter+1;

    for i=1:N

        for k=5:L

            y(k)=X(:,i)'*[y(k-1) y(k-2)*y(k-3) y(k-4) u(k-1)^3 y(k-2)^2 u(k-3)]';

        end

        J=1/(1+(y-y0)*(y-y0)');

        if J>fpbest(i)

            fpbest(i)=J;

            Xpbest(:,i)=X(:,i);

        end 

    end

    [fitnessmax,index]=max(fpbest);

    if fitnessmax>fgbest

        fgbest=fitnessmax;

        Xgbest=X(:,index);

    end

    for i=1:N

        r1=rand; 

        r2=rand;

        fai1=C1*r1;

        fai2=C2*r2;  

        % 速度更新

        V(:,i)=w(iter)*V(:,i)+fai1*(Xpbest(:,i)-X(:,i))+fai2*(Xgbest(:,1)-X(:,i));

        % 若速度超過限定值,則讓其等于邊界值

        index=find(abs(V(:,i))>Vmax);

        if(any(index))

            V(index,i)=V(index,i)./abs(V(index,i)).*Vmax;

        end

        % 位置更新

        X(:,i)=X(:,i)+V(:,i);

    end

%     if (1-fgbest)<err       

%         return  

%     end

    fgbest_fig(iter)=fgbest;

    Xgbest_fig(:,iter)=Xgbest;

end

%--------------------------------------------------------------------figure

% hold on

figure

plot(1:iter_max,fgbest_fig,'-.');

figure

plot(1:iter_max,Xgbest_fig);

%----------------------------------------------------------------------disp

disp(Xgbest)