系统辨识与参数估计大作业第一题递推最小二乘估计参数考虑如上图所示的仿真对象，选择模型结构为：)()2() 1()2() 1()(2121kvkubkubkzakzakz，其中)(kv是服从)1 ,0(N正态分布的不相关随机噪声；输入信号)(ku采用 4 阶逆 M 序列，特征多项式取41)(sssF，幅度为1，循环周期为bitNp62；控制值，使数据的噪信比分别为10％， 73％， 100％三种情况加权因子1)(k；数据长度L=500;初始条件取001. 0)0(?,10)0(6θIP，（1）利用递推最小二乘算法估计参数，（2）利用模型阶次辨识方法（AIC 准则），确定模型的阶次3）估计噪声)(kv的方差和模型静态增益K（4）作出参数估计值随时间的变化图答：设过程的输入输出关系可以描述成( )( )( )Tz khkn k( )z k是输出量，( )h k是可观测的数据向量，n (k)是均值为0 的随机噪声)(kz)(ku)(kv21217.05.115.00.1zzzz217.05.111 zz( )(1),(2),(1), (2)Th kz kz ku ku k1212,,,Ta ab b选取的模型为结构是1212( )(1)(2)(1)(2)z ka z ka z kbu kbu k12121.5,0.7,1.0,0.5aabb加权最小二乘参数估计递推算法RWLS 的公式如下，11( )(1) ( )( ) (1) ( )( )( )(1)( )( )( ) (1)( )( )( )(1)TTTK kp kh khk p kh kkkkK kz khkkp kIK k hkp k为了把 p(k)的对称性，可以把p(k)写成1( )(1)( )( )( ) (1) ( )( )TTp kp kK k Kkhk p kh kk如果把( )k设成 1 的时候，加权最小二乘法就退化成最小二乘法。
用 AIC 准则定阶法来定阶，所用公式nnnnZHV(1), (2), (3),...,( )TnZzzzz L1212,,...,,,...aTnnnaaab bb(0)( 1)...(1)(0)( 1)...(1)(1)(0)...(2)(1)(0)...(2).........................(1)(2)...()(1)(2)...()nzzznuuunzzznuuunHz Lz Lz Lnu Lu Lu Ln其中模型参数n和 噪声( )V k方差的极大似然估计值为ML，2v12()1() ()MLMLMLTT nnnnT vnnnnHHHZZHZHLAIC 的定阶公式写成2 ( )log4vAIC nLn取1,2,3,4;n分别计算( )AIC n， 找到使( )AIC n最小的那个n作为模型的阶次 一般说来，这样得到的模型阶次都能比较接近实际过程的真实阶次信噪比为10%时：参数a1 a2 b1 b2 噪声方差静态增益模型阶次真值-1.5 0.7 1 0.5 1 估计值-1.519 0.72259 1.0314 0.50923 1.0951 7.5661 2 信噪比为73%时：参数a1 a2 b1 b2 噪声方差静态增益模型阶次真值-1.5 0.7 1 0.5 1 估计值-1.519 0.72259 1.0314 0.50923 1.0951 7.5661 2 信噪比为100%时：参数a1 a2 b1 b2 噪声方差静态增益模型阶次真值-1.5 0.7 1 0.5 1 估计值-1.519 0.72259 1.0314 0.50923 估计值7.5661 2 源程序：%function [a1 a2 b1 b2 na nb fangcha Kk]=rwls(L,syn,Np) %na,nb为模型阶次， fangcha 为噪声方差， Kk 为静态增益a1=0;a2=0;b1=0;b2=0;na=0;nb=0;fangcha=0;Kk=0; L=500;Np=62;syn=1; x(1:4)=[1 0 1 0]; for i=1:Np temp=xor(x(1),x(4)); M(i)=x(4); for j=4:-1:2 x(j)=x(j-1); end x(1)=temp; end S=ones(1,Np);%先产生一个全是1 的序列 S if mod(Np,2)==0% 判断 Np 是奇数还是偶数p=Np/2; else p=(Np-1)/2; end for j=1:p S(2*j)=0;% 将 S 序列的偶数位值均置为0，从而使 S 序列是 0 或 1 的方波序列end IM=xor(M,S); % 使用 M 序列与方波序列S 复合生成逆M 序列 IM u=IM*2-1; for i=(Np+1):L u(i)=u(i-Np); end randn('seed',2); v=randn(1,L); syms c; y(1)=0; y(2)=u(1); e(1)=c*v(1); e(2)=1.5*e(1)+c*v(2); for i=3:L y(i)=1.5*y(i-1)-0.7*y(i-2)+u(i-1)+0.5*u(i-2); e(i)=1.5*e(i-1)-0.7*e(i-2)+c*v(i); end m=sum(e.^2); n=sum(y.^2); %c=solve('m/n-syn*syn','c'); c=solve('m/n-syn*syn'); c1=abs(double(c(1))); z(1)=c1*v(1); z(2)=u(1)+1.5*z(1)+c1*v(2); for i=3:L z(i)=1.5*z(i-1)-0.7*z(i-2)+u(i-1)+0.5*u(i-2)+c1*v(i); end cta=zeros(4,L); cta(:,1)=[0.001 0.001 0.001 0.001]'; P=diag([10^6 10^6 10^6 10^6]); %k=2 时h=[-z(1) 0 u(1) 0]'; K= P*h*inv(h'*P*h+1); P=P-K*K'*(h'*P*h+1); cta(:,2)=cta(:,1)+K*(z(2)-h'*cta(:,1)); for k=3:L h=[-z(k-1) -z(k-2) u(k-1) u(k-2)]'; K=P*h*inv(h'*P*h+1); P=P-K*K'*(h'*P*h+1); cta(:,k)=cta(:,(k-1))+K*(z(k)-h'*cta(:,(k-1))); end %以上为参数估计值z=z'; for na=1:4 for nb=1:4 A=zeros(L,na); B=zeros(L,nb); for i=1:L for j=1:na if i>j A(i,j)=z(i-j); end end end for i=1:L for j=1:nb if i>j B(i,j)=u(i-j); end end end H=[A B]; cta1=inv(H'*H)*H'*z;%cta1为模型参数极大似然估计值cgm(na,nb)=(z-H*cta1)'*(z-H*cta1)/L;%cgm为噪声方差极大似然估计值AIC(na,nb)=L*log(cgm(na,nb))+2*(na+nb); end end [na nb]=find(AIC==min(min(AIC))); fangcha=cgm(na,nb)/(c1^2); a1=cta(1,500);a2=cta(2,500);b1=cta(3,500);b2=cta(4,500); Kk=(cta(3,L)+cta(4,L))/(1+cta(1,L)+cta(2,L)); m=1:L; plot(m,cta(1,:),'b-',m,cta(2,:),'k-',m,cta(3,:),'y-',m,cta(4,:),'r-') 第二大题卡尔曼滤波一个系统模型为 )()()1(...1 ,0),()()() 1(22211kwkxkxkkwkxkxkx同时有下列条件：（1）初始条件已知且有。
T] 0, 0[)0(x（2）)(kw是一个标量零均值白高斯序列，且自相关函数已知为jkkwjwE)]()([，另外，我们有下列观测模型，即 ) 1() 1()1(...1 , 0),1() 1() 1(222111 kvkxkzkkvkxkz，且有下列条件：（3）) 1(1kv和) 1(2kv是独立 的零均值白高斯序列，且有jkkvjvE)]()([11，jkkvjvE2)]()([22，...2 ,1 ,0k（4）对于所有的j和k，)(kw与观测噪声过程) 1(1kv和) 1(2kv是不相关的， 即0)]()([1kvjwE，0)]()([2kvjwE我 们 希 望 得 到 由 观 测 矢 量Tkzkzk)]1(),1([) 1(21z估 计 状 态 矢 量Tkxkxk)]1(),1([) 1(21x的卡尔曼滤波的公式表示，并求解以下问题：（a）求 出 卡 尔 曼 增 益 矩 阵 ， 并 得 出 最 优 估 计)1(kx和 观 测 矢 量)1(),...,2(),1 (kzzz之间的递推关系b）用模拟数据确定状态矢量)(kx的估计值)|(?kkx，10...2, 1 ,0k，并画出当10...2, 1 ,0k时，)|(?1kkx，)|(?2kkx的图（c）通常， 状态矢量的真实值是得不到的，但是为了用作图来说明问题，表1 和表 2给出了状态矢量元素的真实值。
对于10...2, 1 ,0k，在同一幅图中画出真实值和在（ b）中确定的)(1kx的估计值对)(2kx重复这一过程当k从 1 变到 10时，对每个元素2, 1i，计算并画出各自的误差图，即)|(?)(kkxkxiid）当k从 1 变到10 时，通过用由卡尔曼滤波器决定的状态误差协方差矩阵画出)]|([2 1kkE和)]|([2 2kkE，而2, 1),|(?)()|(ikkxkxkkiii（e）讨论一下（ C）中计算的误差与（d）中方差之间的关系表 1 观测值时间下标k观测值)(1kz观测值)(2kz1 3． 29691969 2．10134294 2 3． 38736515 0．47540797 3 7． 02830641 3．17688898 4 9． 71212521 2．49811140 5 11．42018315 2．91992424 6 15．97980583 6．17307616 7 22．06934285 5．42519274 8 28．30212781 3．05365741 9 30．44683831 5．98051141 10 38．75875595 4．51016361 表 2 由模拟得到的实际状态值时间下标k实际状态值)(1kx实际状态值)(2kx0 0.00000000 0.00000000 1 1.65428714 1.65428714 2 3.50300702 1.84871988 3 5.99785292 2.47552222 4 9.15040740 3.17187816 5 12.50873910 3.35833170 6 16.92192594 4.41318684 7 21.34483352 4.42290758 8 25.89335144 4.54851792 9 31.54135330 5.64800186 10 36.93605670 5.39447034 答：（a）卡尔曼增益矩阵：kkkTT1HPC(CPCR)’’估计值与观测值之间的递归关系为：kk-1k-1kkkkkXAXH(YCAX)（b）状态矢量估计值的计算框图：（c）)(1kkx和)(2kkx的图：（d）真实值与估计值的比较图：+ 1kH+ 1Z1kA1kC1ky1?kx1?kx1' ?kx各自的误差图：（ e） 通 过 用 卡尔 曼 滤 波 器 的状 态 误 差。