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'); %