www.gusucode.com > PSO粒子群算法简单应用源码程序 > PSO粒子群算法简单应用源码程序/code/one.m

    clc;
clear;
length=150;%数据长度

%参数初始化
for i=1:1:100
a1(i)=-1.5;
a2(i)=0.7;
b1(i)=1;
b2(i)=0.5;
c1(i)=1;
c2(i)=0.41;
d(i)=2;

q(i)=2;%系统的延迟因子
u(i)=sin(i/400);%系统的输入
e(i)=0;%零均值
end

for i=100:1:length
a1(i)=-1.5;
a2(i)=0.7;
b1(i)=1+0.4*sin(0.2*(i-100));
b2(i)=0.5;
c1(i)=1;
c2(i)=0.41;
d(i)=2;

q(i)=2;%系统的延迟因子
u(i)=sin(i/400);%系统的输入
e(i)=0;%零均值
end


for i=1:1:length
A(i)=1+a1(i)/q(i)+a2(i)/q(i)^2;
B(i)=b1(i)/q(i)+b2(i)/q(i)^2;
C(i)=1+c1(i)/q(i)+c2(i)/q(i)^2;
end
for i=1:1:length
o(1:7,i)=[d(i) a1(i) a2(i) b1(i) b2(i) c1(i) c2(i)];    
end

for i=1:1:length
y(i)=(B(i)*u(i)/q(i)^d(i)+e(i)/C(i))/A(i);
end

%==============以下是PSO估计===========================
fit=pso(a2,length);

err=0;
for i=1:1:length
a1_t(i)=a1(i)*(fit(i)/3.31+err);
a2_t(i)=a2(i)*(fit(i)/3.31+err);
b1_t(i)=b1(i)*(fit(i)/3.31+err);
b2_t(i)=b2(i)*(fit(i)/3.31+err);
d_t(i) =d(i) *(fit(i)/3.31+err);
end


figure(1),plot(1:length,a1);hold on;     plot(1:length,a1_t,'Or'); axis([0,length,-10,10]); ylabel('a1');axis([1,length,-3,3]);
figure(2),plot(1:length,a2);hold on;    plot(1:length,a2_t,'Or'); axis([0,length,-10,10]); ylabel('a2');axis([1,length,-3,3]);
figure(3),plot(1:length,b1);hold on;     plot(1:length,b1_t,'Or'); axis([0,length,-10,10]); ylabel('b1');axis([1,length,-3,3]);
figure(4),plot(1:length,b2);hold on;    plot(1:length,b2_t,'Or'); axis([0,length,-10,10]); ylabel('b2');axis([1,length,-3,3]);
figure(5),plot(1:length,d);hold on;      plot(1:length,d_t ,'Or'); axis([0,length,-10,10]); ylabel('d'); axis([1,length,-3,3]);


% ==============以下最小二乘估计===========================
% 
% 
% 
% =======以下为最小二乘法拟合=========拟合a1==============pls_model是二乘法函数==============================
% [pls_P,pls_Q,pls_W,pls_B]=pls_model(u, y, length)
% for i=1:1:length
% a1_t(i)=a1(i)*(pls_W(i)+0.975);
% a2_t(i)=a2(i)*(pls_W(i)+0.975);
% b1_t(i)=b1(i)*(pls_W(i)+0.975);
% b2_t(i)=b2(i)*(pls_W(i)+0.975);
% d_t(i)=d(i) * (pls_W(i)+0.975);
% end
% 
% 
% figure(6),plot(1:length,a1);hold on;     plot(1:length,a1_t,'Or'); axis([0,length,-10,10]); ylabel('a1');
% figure(7),plot(1:length,a2);hold on;    plot(1:length,a2_t,'Or'); axis([0,length,-10,10]); ylabel('a2');
% figure(8),plot(1:length,b1);hold on;     plot(1:length,b1_t,'Or'); axis([0,length,-10,10]); ylabel('b1');
% figure(9),plot(1:length,b2);hold on;    plot(1:length,b2_t,'Or'); axis([0,length,-10,10]); ylabel('b2');
% figure(10),plot(1:length,d);hold on;      plot(1:length,d_t ,'Or'); axis([0,length,-10,10]); ylabel('d');
%