1、系统辨识实验报告姓 名：*学 号：*指导教员：*2011-6-25一、实验目的1. 通过实验，掌握递推最小二乘参数辨识方法。2. 通过实验，掌握残差的平方和 F 检验模型阶次辨识方法。二、实验步骤1. 掌握最小二乘递推算法和 F 检验模型阶次辨识的基本原理。2. 设计实验方案。3. 编制实验程序。4. 调试程序，研究实验问题，记录数据。5. 分析实验结果，完成实验报告。三、实验内容1、实验对象设定实验所用的仿真模型如下： )(2(5.0)1()2(7.0)1(5.)( kvukzkz 其中，u(k) 和 y(k)分别为模型的输入和输出变量；v(k)为零均值、方差为1、服从正态分布的白噪声； 为
2、噪声的标准差（实验时，可取0.0、0.1、0.5、1.0） ；输入变量 u(k)采用 M 序列，幅度取 1.0。辨识模型的形式取 )()11 keuzBzA为方便起见，取 ，即nba nzbzzBa 21)(2、递推最小二乘(RLS)参数辨识算法设计当观测数据长度为 N 时，系统可以写为 NNY递推最小二乘法辨识公式可以写为217.05.z21zy(k)u(k)v(k)+y(k)e(k) NTNTNN PPKy1111 )(其中，初值可以按如下形式给出 00)(TNNY同时计算如下性能指标：1) 损失函数的递推计算 2211()()()kkTii iJzz211()()TNTNJkJkP2)
3、噪声标准差的估计 ()dimJ3) 模型静态增益估计 1niibKa根据上述计算公式可以编制 RLS 的算法步骤如下：step1： 确定被辨识系统的模型结构，以及多项式 和 的阶次；)(1zA)(1Bstep2： 设定递推参数的初值 ， ；)0(Pstep3： 采样获取新的观测数据 和 ，并且组成观测数据向量 ky)(u空间 ；)1(kstep4： 用 RLS 法计算当前参数递推估计值 ；)(kstep5： 采样次数 加 1，然后转回到第三步继续循环。由算法可以编制程序见附录程序 1。3、F-TEST 模型定阶算法设计如果模型为 )()()(111 kzckuzbkyza则残差为 niiezz
4、ze111 )()()()(NnkeJ12)(由于 随着 的增加而减小，在阶数 的增大过程中，我们对那个使 显nJ nJ著减小的阶 感兴趣。为此引入准则1i ()1)2(, (,2)JnNnt FNn其中， 为相应阶次下的损失函数值， 为所用的数据长度， 为模型)J L的估计阶次。若 ,拒绝 ,若 ,接受 ，其中tnt1,( 0:nHtNt1( 001:nH为风险水平 下的阀值。这时模型的阶次估计值可取 。t n同时计算如下性能指标：1) 参数估计平方相对偏差 21 ,niiii2) 参数估计平方根偏差 212,niiii3) 静态增益估计相对偏差 11,Knni ii iKbbaa 可以编制
5、程序见附录程序 2。4、实验结果及分析(1) RLS 参数辨识仿真白噪声 v(k)为零均值、方差为 1、服从正态分布。在程序中用服从正态分布的随机数产生。噪声的标准差 在实验时，分别取 0.0、0.1、0.5、1.0；输入变量 u(k)采用 M 序列,幅度取 1.0，通过异或 xor 移位寄存器循环产生；数据长度 N 可取 100、300、500。0 100 200 30000.20.40.60.81 M估估0 100 200 300-4-2024 估估估估估图 1 输入的 M 序列 图 2 高斯白噪声被辨识参数的初始值 应该是已给充分小的向量；初始状态 应该是)( )(P一个充分大的初始向量
6、 IP2)0(a下图给出了 时的仿真曲线：1.00 100 200 300-3-2-10123Parameter identification with recursive least squares0 100 200 300-3000-2000-10000100020003000 identification precision图 3 参数辨识结果 图 4 参数误差从左图可以大略看出，仿真参数值基本与给定模型一致，而由右图知参数误差非常小。为进一步分析结论，分别取 为 0.0、0.1、0.5、1.0 时，记录参数辨识结果和误差如下表所示：表 1 辨识参数值与误差0 0.1 0.5 11a-1
7、.5 0 -1.4813 9.8e-006 -1.5111 -0.0019 -1.5269 6.9e-00420.7 0 0.6831 1.7e-005 0.6908 -0.0031 0.7364 0.0010b1 0 0.9967 -5.5e-005 0.9461 0.0039 1.0842 0.00670.5 0 0.5206 -1.3e-004 0.4202 0.0148 0.5982 0.0104注：在数据单元内，左边为仿真参数值，右 边为参数误差。从以上表格可以看出，辨识精度较高，基本上满足要求；同时，随着噪声方差的增加，辨识精度也基本上长呈降低的趋势。同时，计算不同噪声方差下的性能
8、指标如下：表 2 参数辨识性能指标0 0.1 0.5 1噪声标准差 0.6809 1.5731 4.6072 7.6446模型静态增益 K0 7.4886 6.9125 6.4417损失函数 J44.5082 237.5618 2.0377e+003 5.6103e+003由上表可以看出，随着噪声的加大，系统精度也降低。根据图 3、图 4 和表 1、表 2，也不难发现递推最小二乘法算法的不足： 在上述递推算法过程中，所有先前数据的影响一致在起作用； 随着递推运算的逐步深入，参数估计将收敛到某固定值，因此，适合于定常未知参数的估计； 在时变参数的情况下，参数时变信息更多地蕴藏在新观测数据中，而与
9、先前观测数据的关系将逐渐减弱，因此，它不适合时变参数的估计。(2) F-TEST 模型定阶仿真根据程序 2，画出模型阶次曲线如下图所示：1 1.5 2 2.5 3 3.5 40.511.522.5 估估估估估估估估估估估估估估估估估估估估图 5 曲线图nJ由图 5 可以看出，当 为 2 时 最后一次出现陡峭的下降， 再大，则nnJn只有微小的变化。nJ由程序2可得 的值，根据如下公式计算的 的值，利用各阶残差的平方和Jt进行F检验，当 时， F的阀值为3.09，构建 F检验辨识系统阶次表如表305.所示： 2)1()1,( nNnJnt表3 F检验表n1 2 3 4J2.4406 0.9568
10、 0.9542 0.9535t777.72 1.36 0.337因此可以选定模型阶次为 2。n附录：程序清单程序 1：rls 参数辨识 clc;clear;L=100; %M序列异或发生器x1=1;x2=1;x3=1;x4=0; for k=1:L; u(k)=xor(x3,x4); x4=x3;x3=x2;x2=x1;x1=u(k);end %产生高斯白噪声v=randn(1,L);q=0;%初始值设定c0=0.001 0.001 0.001 0.001; p0=106*eye(4,4);E=0.000000005; E1=-0.000000005;c=c0,zeros(4,L-1);e=z
11、eros(4,L);y=zeros(1,L);ji=zeros(1,L);tt=0;jie=0;pp0=p0;%最小二乘递推确定模型参数for k=3:L;y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-1)+0.5*u(k-2)+q*v(k); h1=-y(k-1),-y(k-2),u(k-1),u(k-2);x=h1*p0*h1+1;x1=inv(x);k1=p0*h1*x1;d1=y(k)-h1*c0;c1=c0+k1*d1;e1=c1-c0;e2=e1./c0;e(:,k)=e2;c0=c1; c(:,k)=c1;ji(k)=d12;%=损失函数的计算=%for m=1:
12、k-1tt=tt+ji(m);endjie=tt+(y(k)-h1*c(:,k-2)2/(h1*pp0*h1+1);%=%p1=p0-k1*k1*h1*p0*h1+1;pp0=p0;p0=p1; if e2E1 break;end endjiea1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);A1=a1(L) A2=a2(L)B1=b1(L) B2=b2(L)%模型静态增益kjian=(B1+B2)/(1+A1+A2)%噪声标准差估计lanmuda=sqrt(jie/(L-4)
13、subplot(2,1,1);i=1:L;plot(i,a1,r,i,a2,:,i,b1,g,i,b2,:)title(Parameter identification with recursive least squares method)subplot(2,1,2);i=1:L; plot(i,ea1,r,i,ea2,g,i,eb1,b,i,eb2,r)title(identification precision);程序 2：f-test 模型定阶%=close allclearclc%=产生M序列作为输入 =x=0 1 0 1 1 0 1 1 1; %初始值n=1003; %n为脉冲数目
14、M=; %存放M 序列for i=1:ntemp=xor(x(4),x(9);M(i)=x(9);for j=9:-1:2x(j)=x(j-1);endx(1)=temp;end%产生高斯白噪声v=randn(1,1004);z=;z(1)=-1;z(2)=0;L=1000;for i=3:L+4z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i);end%模型阶次n=1for i=1:LH1(i,1)=z(i);H1(i,2)=M(i);endestimate=inv(H1*H1)*H1*z(2:1001);e=z(2:1001)-H1*esti
15、mate;V1=e*e/L%模型阶次n=2for i=1:LH2(i,1)=z(i+1);H2(i,2)=z(i);H2(i,3)=M(i+1);H2(i,4)=M(i);endestimate=inv(H2*H2)*H2*z(3:1002);e=z(3:1002)-H2*estimate;V2=e*e/L% n=3for i=1:LH3(i,1)=z(i+2);H3(i,2)=z(i+1);H3(i,3)=z(i);H3(i,4)=M(i+2);H3(i,5)=M(i+1);H3(i,6)=M(i);endestimate=inv(H3*H3)*H3*z(4:1003);e=z(4:1003)-H3*estimate;V3=e*e/L%模型的阶次n=4for i=1:LH4(i,1)=z(i+3);H4(i,2)=z(i+2);H4(i,3)=z(i+1);H4(i,4)=z(i);H4(i,5)=M(i+3);H4(i,6)=M(i+2);H4(i,7)=M(i+1);H4(i,8)=M(i);endestimate=inv(H4*H4)*H4*z(5:1004);e=z(5:1004)-H4*estimate;V4=e*e/Lplot(1:4,V1 V2 V3 V4)title(利用残差的方差估计模型的阶次 )xlabel(阶次 )ylabel(残差方差 )
