标准PSO辨识NARMAX模型源码程序

- 日理万妓 2022-09-04 12:41 39阅读 0赞

% 用标准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)

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2Njc3NzMjI_size_16_color_FFFFFF_t_70][]

B-55

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2Njc3NzMjI_size_16_color_FFFFFF_t_70]: /images/20220829/72fd8ff209314a6aab3b888f5f348c08.png