psojb.m

利用基本的PSO算法对传感器进行建模

clear all;
clc;

fs=2500;N=2500;
t=0:1/fs:(N-1)/fs;
global y yd

% 传感器的二阶模型的阶跃响应 2500点
num_cgq=[7318];
den_cgq=[1,127.1,7318];
sys_cgq=tf(num_cgq,den_cgq);
[y]=step(sys_cgq,t);
%是补偿器输入

%等效系统的阶跃响应 2500点
num_lx=[18225 ]  ;den_lx=[ 1, 190.89,18225 ]; 
sys_lx=tf(num_lx,den_lx);
zpk(sys_lx);
[yd]=step(sys_lx,t);
%size(y_lx)
%yd(1:2500)是参考模型输出

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

% 随机初始化种群的位置和速度
x=randn(popsize,dimen);

%%x(1,:)为局部最优
v=(rand(popsize,dimen)-0.5)*2;              % 随机产生速度范围[-1,1]
vmax=ones(1,dimen);
% 先计算每个粒子的适应度，并初始化个体最优Pbest和群体最优Pgbest
for i=1:popsize
    pbest_fit(i)=fitness(x(i,:),dimen);
    xpbest(i,:)=x(i,:);                %  初始化个体最优Pbest
end
                      
[gbest_fit index]=min(pbest_fit(i));      % 初始化群体最优Pgbest的初始值
temp_gbest_fit=gbest_fit;
xgbest=x(index,:);

%进入主循环，直到迭代结束
for diedai=1:L
    disp(['No.Loop : ' num2str(diedai) ]) ;
       for i=1:popsize
        for j=1:dimen
            v(i,j)=w*v(i,j)+c1*rand(1)*(xpbest(i,j)-x(i,j))+c2*rand(1)*(xgbest(j)-x(i,j));
            if abs(v(i,j))>vmax(j)             % 速度限制范围
                v(i,j)=sign(v(i,j))*vmax(j);
            end
            x(i,j)=x(i,j)+v(i,j);
        end
        temp_x_fit=fitness(x(i,:),dimen);
        if temp_x_fit<pbest_fit(i)
            pbest_fit(i)=temp_x_fit;
            xpbest(i,:)=x(i,:);    %  第i个粒子个体最优更新
        end
        if pbest_fit(i)<temp_gbest_fit
            temp_gbest_fit=pbest_fit(i);
            xgbest=xpbest(i,:);    % 用i粒子的个体最优来更新群体最优
        end
    end
    gbest_fit(diedai)=temp_gbest_fit;      % 记录第*次迭代得到的全局最好适应度函数
end

disp('全局最优位置为：')
xgbest'
disp('全局最优适应度函数值：')
gbest_fit(diedai)

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

%带入实验数据，看滤波后的效果
yc(1:3)=0;
for k=4:2500
    yc(k)=xgbest(1)*y_data(k)+xgbest(2)*y_data(k-1)+xgbest(3)*y_data(k-2)+xgbest(4)*y_data(k-3)-xgbest(5)*yc(k-1)-xgbest(6)*yc(k-2)-xgbest(7)*yc(k-3);
end
%[b2,a2]=butter(2,0.01);
%yc=filter(b2,a2,yc);
figure(1)
plot(t,yc,'k',t,y_data,'k')
xyh_plot
axis([0,1,-0.5,2.5]);
title('flux=50g/s')
annotation1 = annotation(figure(1),'line',[0.5335 0.492],[0.8127 0.7745]);
annotation1 = annotation(figure(1),'line',[0.5716 0.598],[0.7314 0.6644]);
annotation1 = annotation(...
  figure(1),'textbox',...
  'Position',[0.5902 0.5305 0.3396 0.1489],...
  'LineStyle','none',...
  'FontName','Times New Roman',...
  'FontSize',8,...
  'String',{'1'},...
  'FitHeightToText','on');
annotation1 = annotation(...
  figure(1),'textbox',...
  'Position',[0.4448 0.6426 0.3396 0.1489],...
  'LineStyle','none',...
  'FontName','Times New Roman',...
  'FontSize',8,...
  'String',{'2'},...
  'FitHeightToText','on');

figure(3)
plot(gbest_fit,'k') % 适应值的变化
xlabel('迭代次数(L)','fontname','Times New Roman','Fontsize',8);
ylabel('适应度函数(V)','fontname','Times New Roman','Fontsize',8);
set(gcf,'color',[1,1,1]);
set(gca,'xcolor',[0,0,0],'ycolor',[0,0,0]);
set(gcf,'units','centimeters','position',[5,10,6.8,5.2]);
set(gca,'box','on','fontname','Times New Roman','fontsize',8);