%%%%%%%%%%%以上为数据滤波%%%%%%%% phi(:,1)=-yy(n:n+N-1); phi(:,2)=-yy(n-1:n+N-2); phi(:,3)=uu(n:n+N-1); phi(:,4)=uu(n-1:n+N-2); theta=inv(phi'*phi)*phi'*yy(n+1:n+N); end f theta ELS %增广最小二乘的递推算法 %===========产生均值为0，方差为1 的高斯白噪声========= v(1)=0; v(2)=0; z=UY(:,2); M=UY(:,1); %递推求解 P=100*eye(6); %估计方差 Pstore=zeros(6,800); Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)]; Theta=zeros(6,401); %参数的估计值，存放中间过程估值 Theta(:,1)=[3;3;3;3;3;3]; K=[10;10;10;10;10;10]; for i=3:801 h=[-z(i-1);-z(i-2);M(i-1);M(i-2);v(i-1);v(i-2)]; K=P*h*inv(h'*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2)); v(i)=z(i)-h'*Theta(:,i-2); P=(eye(6)-K*h')*P; Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)]; end %======================================================================= disp('参数a1、a2、b1、b2、d1、d2 估计结果：') Theta(:,800)
RLS
三阶
-0.4381-0.12280.0407 -0.00780.5652 0.5106 0.1160
以下给出RLS中的参数估计过程曲线和误差曲线
程序清单： LS（二阶）： M=UY(:,1); z=UY(:,2); H=zeros(199,5); for i=1:199 H(i,1)=-z(i+1); H(i,2)=-z(i); H(i,3)=M(i+2); H(i,4)=M(i+1); H(i,5)=M(i); end Estimate=inv(H'*H)*H'*(z(3:201)) RLS（二阶）： clc M=UY(:,1); z=UY(:,2);
自动化09-3
宋佳瑛
09051304
系统辨识实验报告 实验一：系统辨识的经典方法 实验目的：掌握系统的数学模型与系统的输入，输出信号之间的关系， 掌握经典辨识的实验测试方法和数据处理方法。熟悉matlab/Simulink 环境。 实验内容： 1.用系统阶跃响应法测试给定系统的数学模型。 在系统没有噪声干扰的条件下通过测试系统的阶跃响应获得系统的一阶 加纯滞后或二阶加纯滞后模型，对模型进行验证。 2.在被辨识的系统加入噪声干扰，重复上述1的实验过程。 1. 没有噪声 搭建对象
实验程序：GLS %%%%%%%%%%%%%%%%广义最小二乘法辨识%%%%%%%%%%% clc n=2; N=799; % y=y1; % u=u1; y=UY(:,2); u=UY(:,1); % y=y3; % u=u3; Y=y(n+1:n+N); phi1=-y(n:n+N-1); phi2=-y(n-1:n+N-2); phi3=u(n:n+N-1); phi4=u(n-1:n+N-2); phi=[phi1 phi2 phi3 phi4]; theta=inv(phi'*phi)*phi'*Y %%%%%%以上为最小二乘计算Theta估计值%%%%%%%%% f0=[10;10]; while 1 e(n+1:n+N,1)=y(n+1:n+N)-phi*theta; omega(:,1)=-e(n:n+N-1); omega(:,2)=-e(n-1:n+N-2); f=inv(omega'*omega)*omega'*e(n+1:n+N,1); if (f-f0)'*(f-f0)<0.0001 break end f0=f; %%%%%%以上为最小二乘计算f估计值%%%%%%%%% for i=n+1:n+N yy(i,1)=[y(i),y(i-1),y(i-2)]*[1,f(1),f(2)]'; uu(i,1)=[u(i),u(i-1),u(i-2)]*[1,f(1),f(2)]'; end yy(1,1)=y(1,1); uu(2,1)=u(2,1);
测试对象流程图 实验结果为：
2、加入噪声干扰 搭建对象
实验结果：
加入噪声干扰之后水箱输出不平稳，有波动。
实验二：相关分析法 搭建对象：
处理程序： for i=1:15 m(i,:)UY(32-i:46-i,1);
end y=UY(31:45,2); gg=ones(15)+eye(15); g=1/(25*16*2)*gg*m*y; plot(g); hold on; stem(g); 实验结果： 相关分析法
P=100*eye(5); %估计方差 Pstore=zeros(5,200); Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5)]'; Theta=zeros(5,200); %参数的估计值，存放中间过程估值 Theta(:,1)=[0;0;0;0;0]; K=[10;10;10;10;10;10;10]; for i=3:201 h=[-z(i-1);-z(i-2);M(i);M(i-1);M(i-2)]; K=P*h*inv(h'*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2)); P=(eye(5)-K*h')*P; Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5)]'; end i=1:200; figure(1) plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)) title('待估参数过渡过程') figure(2) plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,: title('估计方差变化过程') 同理可以写出三阶的LS以及RLS算法，此处略去。 2.在t=250出加入一个阶跃扰动，令扰动值为0.05 利用RLS和带遗忘因子的RLS计算结果如下表。 算法 RLS RLS 遗忘因 A1 A2 A3 B0 B1 B2 B3 子 1 -0.4791-0.09170.0349 -0.00680.5659 0.4882 0.1033 0.98 -0.5498-0.03900.0270 -0.00690.5628 0.4515 0.0777 RLS的参数变化过程曲线和误差曲线如下：
i=1:800; figure(1) plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:),i,T title('待估参数过渡过程') figure(2) plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,: title('估计方差变化过程')
带0.98遗忘因子的参数过度曲线和误差曲线
根据上面两种方法所得到的误差曲线和参数过渡过程曲线，我们可以看
出来利用最小二乘法得到的参数最终趋于稳定，为利用带遗忘因子的最 小二乘算法，曲线参数最终还是有小幅度震荡。由此可以看出两种算法 的一些特点与区别。 最小二乘法：递推算法没获得一次新的观测数据就修正一次参数估计 值，随着时间的推移，便能获得满意的辨识结果。 带遗忘因子的最小二乘法。其本质还是最小二乘法，只不过加强新的数 据提供的信息量，逐渐削弱老的数据，防止数据饱和。 2. 参照index3，设计符合GLS和ELS的对象模型，改写参照程序， 实现相应的算法。 GLS 与ELS的实现 参照index3，设计符合GLS和ELS的对象模型，改写参照程序，实现相应 的算法。 利用GLS算法求出参数如下： a1 a2 b1 b2 -0.7411 0.0811 0.2364 0.1040 即y(k)=0.7411y(k—1)-0.0811y(k-2)+0.2364u(k-1)+0.1040u(k-2); 利用ELS求解，得到如下的结果 参数 a1 a2 b1 b2 d1 d2 ELS -0.7424 0.0819 0.2346 0.1039 0 0 绘制出ELS算法下参数的变化曲线以及误差曲线如下：
第二题 aa=5;NNPP=7;ts=2; RR=ones(7)+eye(7); UU= [UY(15:21,1)';UY(14:20,1)';UY(13:19,1)';UY(12:18,1)';UY(11:17,1)'; UY(10:16,1)';UY(9:15,1)']; YY=[UY(8:14,2)]; GG=(RR*UU*YY+4.4474)/[aa*aa*(NNPP+1)*ts]; plot(0:2:13,GG); hold on stem(0:2:13,GG,'filled'); grid on; title('脉冲响应序列') 求系统的脉冲传递函数 g=[ -0.0036 0.2922 0.0820 0.0489]'; a=[]; for k=1:1:6 a=[a;g(k),g(k+1)]; end G=[a;g(7),g(1)]; G1=[g(3:7,1);g(1);g(2)]; A=G\(-G1); a1=A(2) a2=A(1) D=[1 0;a1 1]; C=[g(1);g(2)]; B=D*C; b1=B(1) b2=B(2) zuixiao 相关二步法 z=UY(:,2); M=UY(:,1);