psojmfm2.m

将动态标定实验的阶跃响应数据进行处理之后

%传感器模型为2阶，用PSO
clear all;
clc;

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

% 阶跃输入给阀门
t1=1171;                        
for i=1:t1
     x1(i)=0;
end
for i=(t1+1):2500
    x1(i)=1;
end
x1=x1';

%流量幅值输入到非线性静态环节
winput=50.0467;
p=[0.00000000158458,-0.00000044392527,0.00004927931552,-0.00288655624523,0.10954571440471,0.94145172599242];
winput=polyval(p,winput)-0.94145172599242;
x=winput*x1;% 阶跃电压输入给阀门
%阀门传递函数
num_fm=[1]  ;
den_fm=[0.012,1];
[num_fmls,den_fmls] = c2dm(num_fm,den_fm,1/2500,'method');%转化为阀门离散模型
y_fm=filter(num_fmls,den_fmls,x); %阀门输出阶跃响应，二阶传感器模型的输入x(k)

%********************
%传感器输出的实验数据,即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=filter(b,a,y_data);
%y(1:2500)是补偿器输入

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

% 随机初始化种群的位置和速度

b1=-0.001*rand(popsize,1);
b2=0.005*rand(popsize,1);
a1=-2*rand(popsize,1);
a2=rand(popsize,1);
x=[b1,b2,a1,a2];

v=(rand(popsize,dimen)-0.5)*6;              % 随机产生速度范围[-3,3]
vmax=3*ones(1,dimen);
% 先计算每个粒子的适应度，并初始化个体最优Pbest和群体最优Pgbest
for i=1:popsize
    pbest_fit(i)=fit2(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) ]) ;
    w(diedai)=(wmax-wmin)*(L-diedai)/L+wmin;
    c2=0.8+diedai/L;
    for i=1:popsize
        for j=1:dimen
            v(i,j)=w(diedai)*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=fit2(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)


%带入实验数据，看模型的输出和实际输出是否吻合
yc(1:2)=0;
for k=3:2500
    yc(k)=xgbest(1)*y_fm(k-1)+xgbest(2)*y_fm(k-2)-xgbest(3)*yc(k-1)-xgbest(4)*yc(k-2);
end
figure(1)
plot(t,yc,'k',t,y,'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(2)
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);

figure(3)
plot(gbest_fit(50:L),'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);

gbest_fit(1:50);