psobc318.m

利用PSO算法设计传感器的动态补偿器

clear all;
clc;
% 15行7列
popsize=15;  % 种群规模，即粒子数15个,
dimen=7;         % 维数，即要求7个参数    x=[b0,b1,b2,b3,a1,a2,a3]
L=100;      % 最大迭代次数,         区别大小写
error=10^(-6); %迭代终止的精度要求
wmax=0.9;      % 最大惯性权重
wmin=0.4;    % 最小惯性权重
c1=2;c2=1.6;   % 学习因子
w=0.7298;

% 随机初始化种群的位置和速度
for i=1:popsize
    for j=1:dimen
        x(i,j)=randn;
        v(i,j)=randn;
    end
end

% 先计算每个粒子的适应度，并初始化个体最优Pbest和群体最优Pgbest
for i=1:popsize
    pbest_fit(i)=fitness(x(i,:),dimen);
    pbest(i,:)=x(i,:);                %个体最优Pbest
end
pgbest=x(1,:);                       % 群体最优Pgbest的初始值
for i=2:popsize
    if fitness(x(i,:),dimen)<fitness(pgbest,dimen)
        pgbest=x(i,:);               % 群体最优Pgbest
    end
end
temp_gbest_fit=fitness(pgbest,dimen);
%进入主循环，直到迭代结束
for diedai=1:L
    disp(['No.Loop : ' num2str(diedai) ]) ;
    for i=1:popsize
        v(i,:)=w*v(i,:)+c1*rand*(pbest(i,:)-x(i,:))+c2*rand*(pgbest-x(i,:));
        x(i,:)=x(i,:)+v(i,:);
        temp_x_fit=fitness(x(i,:),dimen);
        if temp_x_fit<pbest_fit(i)
            pbest_fit(i)=temp_x_fit;
            pbest(i,:)=x(i,:);    %第i个粒子个体最优更新
        end
        if pbest_fit(i)<temp_gbest_fit
            temp_gbest_fit=pbest_fit(i);
            pgbest=pbest(i,:);    % 用i粒子的个体最优来更新群体最优
        end
    end
    pgest_fit(diedai)=temp_gbest_fit;      %第*次迭代得到的全局最好适应度函数
end

disp('全局最优位置为：')
pgbest'
disp('全局最优适应度函数极值为：')
result=fitness(pgbest,dimen)

%传感器输出的实验数据,即y(1:2500)是补偿器输入
fp=fopen('liu50.txt','r');
out=fscanf(fp,'%f',2500);
fclose(fp);
for i=1:2500
    y(i)=out(i)-0.94404;
end
y=y';
[b,a]=butter(2,0.1);
y=filter(b,a,y);
%y(1:2500)是补偿器输入

%理想连续传递函数
num_lx=[90000 ]  ;
den_lx=[ 1,424.2,90000 ];
sys_lx=tf(num_lx,den_lx);
zpk(sys_lx);
[num_lsan,den_lsan]=c2dm(num_lx,den_lx,1/2500,'method');%
t1=1164;
for i=1:t1
    u_in(i)=0;
end
for i=(t1+1):2500
    u_in(i)=2.14;
end
u_in=u_in';
yd=filter(num_lsan,den_lsan,u_in);
%yd(1:2500)是参考模型输出

yc(1:3)=0;
for k=4:2500
    yc(k)=pgbest(1)*y(k)+pgbest(2)*y(k-1)+pgbest(3)*y(k-2)+pgbest(4)*y(k-3)-pgbest(5)*yc(k-1)-pgbest(6)*yc(k-2)-pgbest(7)*yc(k-3);
end
fs=2500;N=2500;
t=0:1/fs:(N-1)/fs;
figure(1)
plot(t,yc,'m',t,yd,'r',t,y,'k')
xyh_plot