在Hammerstein非线性模型中，基于PSO的参数辨识系统

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