Hammerstein非线性模型的基于PSO的参数辨识系统的本质就是将参数的辨识问题转换为参数空间优化问题，对整个参数域进行搜索并最终获得最优的参数估计。我们需要的参数辨识模型具体描述如下所示：
将Hammerstein非线性模型进行分离，得到8个不同的模型，逐个对其参数进行识辨。下面从较为简单的NL8开始说明直到NL1，NL这九个结果。最后对整个Hammerstein非线性模型进行识辨仿真。
NL
原非线性模型的坐标图。
PSO最佳适应曲线。
V（t）中七个参数的辨识过程。
Y（t）的三个参数的辨识过程如下所示：
部分核心代码如下：
clc;
close all;
clear;
%对a1，b1，b2进行参数辨识
Len = 1000;
%开始粒子群优化
%开始粒子群优化
%粒子群参数设置
ys = zeros(1,Len);
y0 = zeros(1,Len);
var = 0.001; % 噪声方差
err = 0.01;
iter_max = 200; % 最大迭代次数
N = 50; % 种群规模
C1 = 2; % 加速度常数
C2 = 2;
Xmin = -1; % 解取值范围[Xmin,Xmax]
Xmax = 1;
p = 3; % 粒子维数
w = linspace(0.9,0.4,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 = 1;
V = Vmax*(2*rand(p,N)-1);
e = idinput(Len,'rgs',[0 1],[-var var]); % 高斯白噪声，均值为0，方差为var
%定义输入输出
load data1.mat
clear u alpha v f1 f2 f3 h1 h2 h1pr h2pr
load result.mat
clear Xgbest_fig p_best pa_best pb_best Z1_best Z2_best e1_best e2_best Xgbest_fig
for t=1:Len
y0(t) = y(t) + e(t);
end
iter=0;
while iter<iter_max
iter=iter+1;%进行迭代
iter
for i=1:N
for t = 3:Len
a1 = X(1,i);
a2 = X(2,i);
b1 = X(3,i);
%y(t)
ys(t) = -a1*ys(t-1) - a2*ys(t-2) + vs(t-1) + b1*vs(t-2);
end
J=1/(1+(ys-y0)*(ys-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
fgbest_fig(iter) = fgbest;%种群最佳适应度值
Xgbest_fig(:,iter) = Xgbest;%种群最佳位置
end
figure;
plot(1:iter_max,fgbest_fig,'r-*');
%参数优化过程曲线
figure
plot(1:iter_max,Xgbest_fig(1,:),'r-.','LineWidth',2);
hold on
plot(1:iter_max,Xgbest_fig(2,:),'b-.','LineWidth',2);
hold on
plot(1:iter_max,Xgbest_fig(3,:),'m-.','LineWidth',2);
hold on
legend('a1 PSO','a2 PSO','b1 PSO');
plot(1:iter_max,a1,'k');
hold on
plot(1:iter_max,a2,'k');
hold on
plot(1:iter_max,b1,'k');
hold off
disp(Xgbest)
a1_best = Xgbest(1);
a2_best = Xgbest(2);
b1_best = Xgbest(3);
figure;
subplot(5,2,8);
plot(1:iter_max,Xgbest_fig(1,:),'r-.','LineWidth',2);
hold on
plot(1:iter_max,a1,'k');
hold off
xlabel('迭代次数');
ylabel('a1')
subplot(5,2,9);
plot(1:iter_max,Xgbest_fig(2,:),'r-.','LineWidth',2);
hold on
plot(1:iter_max,a2,'k');
hold off
xlabel('迭代次数');
ylabel('a2')
subplot(5,2,10);
plot(1:iter_max,Xgbest_fig(3,:),'r-.','LineWidth',2);
hold on
plot(1:iter_max,b1,'k');
hold off
xlabel('迭代次数');
ylabel('b1')
load result.mat
subplot(5,2,1);
plot(1:iter_max,Xgbest_fig(1,:),'r-.','LineWidth',2);
hold on
p = 2;
plot(1:iter_max,p,'k');
hold off
xlabel('迭代次数');
ylabel('p');
subplot(5,2,2);
plot(1:iter_max,Xgbest_fig(2,:),'r-.','LineWidth',2);
hold on
pa = 0.5;
plot(1:iter_max,pa,'k');
hold off
xlabel('迭代次数');
ylabel('pa')
subplot(5,2,3);
plot(1:iter_max,Xgbest_fig(3,:),'r-.','LineWidth',2);
hold on
pb = 1;
plot(1:iter_max,pb,'k');
hold off
xlabel('迭代次数');
ylabel('pb')
subplot(5,2,4);
plot(1:iter_max,Xgbest_fig(4,:),'r-.','LineWidth',2);
hold on
Z1 = 1.2;
plot(1:iter_max,Z1,'k');
hold off
xlabel('迭代次数');
ylabel('Z1')
subplot(5,2,5);
plot(1:iter_max,Xgbest_fig(5,:),'r-.','LineWidth',2);
hold on
Z2 = -1;
plot(1:iter_max,Z2,'k');
hold off
xlabel('迭代次数');
ylabel('Z2')
subplot(5,2,6);
plot(1:iter_max,Xgbest_fig(6,:),'r-.','LineWidth',2);
hold on
e1 = 1.5;
plot(1:iter_max,e1,'k');
hold off
xlabel('迭代次数');
ylabel('e1')
subplot(5,2,7);
plot(1:iter_max,Xgbest_fig(7,:),'r-.','LineWidth',2);
hold on
e2 = -1.2;
plot(1:iter_max,e2,'k');
hold off
xlabel('迭代次数');
ylabel('e2')
A27-01
